Mathematics of Computation 


A Quarterly Journal 


Edited by 


ALAN FLETCHER P. C. HAMMER EUGENE ISAACSON 
Y. L. LUKE DANIEL SHANKS A. H. TAUB 
R. 8. VARGA J. W. WRENCH, JR. D. M. YOUNG 


HARRY POLACHEK, Chairman 


NOS. 69-72 


1960 


Formerly: Mathematical Tables and other Aids to Computation 


Published by the 
National Academy of Sciences—National Research Council 


Washington, D. C. 








/ 
ron 








The SIAM Review 








ARTICLES 


On the role of professional societies in stimulating and guiding of research 

R. F. Drenick 

Optimum allocation of discharge to units in a hydroelectric generating 

IN poor a Lig sshd Ws Wig drink Osi Klose Siok ASS AO B. Bernholtz 
Applications of graphs and boolean matrices to computer programming 

Rosalind B. Marimont 





A note concerning orthogonal polynomials.......... Raiph Hoyt Bacon 
Special block iterations with applications to Laplace and biharmonic dif- 

os raed s cons nace es aes baw 4 .. Herbert B. Keller 
A theorem on determinants...................-.-. ........Huyen Gott 
Ni Bic seda caneeisensdnowasnae ee E. 8S. Wolk 

PROBLEMS 

kn s 6k eRe eee oe Herbert A. Cohen 
ey Catan sc a Kins pais Rabin ew d cw S aaaed Larry Shepp 
a te sh a racist admis lee nacd D. J. Newman 
en ink 5 dk nae duds neu oh + ke acneme W. L. Bade 


A sorting problem 


BOOK REVIEWS 


Approximate Methods of Higher Analysis (Kantorovich and Krylov) 
T. J. Rwlin 
Handbook of Automation, Computation, and Control, Volume 2, Com 
puters and Data Processing (Grabbe, Ramo, and Wooldridge) 
Richard M. Karp 


NEWS AND NOTICES 














January 1960 - Vol. 14, No. 69 


Mathematics... 


FEB -3 1060 
of oa 


Computation 


A journal devoted to advances in numerical analysis, 
the application of computational methods, mathematical tables, 
high-speed calculators and other aids to computation 





Formerly: Mathematical Tables and other Aids to Computation 


Published Quarterly by the 
National Academy of Sciences—National Research Council 





Editorial Committee 
Division of Mathematics 
National Academy of Sciences—National Research Council 
Washington, D. C. 


H. Po.acHex, oi a ae Mathematics Laboratory, David Taylor Model 
Basin, Washington 7, 

C. C. Craia, University er Michigan, Ann Arbor, Michigan 

ALAN Fiercuer, University of Liverpool, Liverpool 3, England 

Evucene Isaacson, New York University, New York 3, New York 

D. ar — Mathematics Laboratory, David Taylor Model Basin, Wash- 
ington 

+e pe * Sirs, Goddard Space Flight Center, National Aeronautics and Space 
Administration, Washington 25, D.C. 

A. H. Tavs, University of Illinois, Urbana, Illinois 

C. B. Tompxins, University of California, Los Angeles, California 

J. W. Wrenca, Jr., Applied Mathematics Laboratory, David Taylor Model Basin, 
Washington 7, D. C. 


Information to Subscribers 


The journal is published quarterly in one volume per year with issues numbered 
serially since Volume I, Number 1. Starting with January, 1959 subscriptions are 
$8.00 per year, single copies $2.50. Back issues are available as follows: 
Volume I (1943-1945), Nos. 10 and 12 only are available; $1.00 per issue. 
Volume II (1946-1947), Nos. 13, 14, 17, 18, 19, and 20 only available; $1.00 
per issue. 
Volume III (1948-1949), Nos. 21-28 available. $4.00 per year (four issues), 
$1.25 per issue. 
Volume IV (1950 through 1958), all issues available; $5.00 per year, $1.50 per 
issue. 


Microcard Edition 


Volumes I—X (1943-1956), Nos. 1-56 are now available on Microcards and may be 
from the Microcard Foundation, Box 2145, Madison 5, Wisconsin, at a 
cost of $20.00 for the complete set. Succeeding volumes are available on request. 


Information to Contributors 


All contributions intended for publication in Mathematics of Computation and all 
books for review should be addressed to H. Polachek, Technical Director, Applied 
Mathematics Laboratory, David Taylor Model Basin, Washington 7, D. C. The 
author may suggest an appropriate editor for his paper. Manuscripts should be 
typewritten double-spaced in the format used by the journal. Authors should sub- 
mit the original and one copy, and should retain one copy. 


Subscriptions, address changes, business communications and payments should be 
sent to: 


THE PRINTING AND PUBLISHING OFFICE 
Tue NationaL ACADEMY OF SCIENCES 
2101 Constitution Avenue 

Washington 25, D. C. 


PUBLISHED BY THE 
NATIONAL ACADEMY OF SCIENCES—NATIONAL RESEARCH COUNCIL 
Mt. Royal and Guilford Avenues, Raltimore 2, Md. 


Second-class postage paid at Baltimore, Md. 








Wy 


Ad 





472 


Bo, 6Y5EG MATHEMATICS OF COMPUTATION 


Beginning with the January 1960 issue, the title of Mathematical Tables and 
other Aids to Computation has been changed to Mathematics of Computation. The 
new name was unanimously adopted by the Editorial Committee in order to re- 
flect the broadened scope of the Journal, which has expanded to meet the need 
in this country for a publication devoted to numerical analysis and computation. 
As stated in the announcement of the October issue, the new title in no way repre- 
sents a diminished interest in mathematical tables, which will continue to be a sub- 
ject of major emphasis, as in the past. It recognizes, however, an increased in- 
terest in other areas in the field of Mathematics of Computation, which have grown 
in importance and in which rapid advances are being made in the present era of 
technological progress. 

The first issue of Mathematical Tables and other Aids to Computation appeared 
in January 1943, more than one year prior to the completion of the first large 
scale computer—the Harvard Automatic Sequence Controlled Calculator. It is 
quite evident, however, that Professor R. C. Archibald, when he founded the jour- 
nal, had to a large extent foreseen the great progress which was to follow in the de- 
velopment of instruments for computation. In the introductory remarks which he 
published in the first issue he established as one of the principal goals of the Jour- 
nal, ‘‘to serve as a clearing-house for information” concerning tools of computation 
which have been “vastly multiplied during the past decade.’”’ He further speaks of 
these tools, or accounts of them, which “are to be found in an enormous interna- 
tional range of book, pamphlet, and periodical publication, not only in the fields of 
Pure Mathematics, Physics, Statistics, Astronomy, and Navigation, but also in 
such fields as Chemistry, Engineering, Geodesy, Geology, Physiology, Economics, 
and Psychology.” 

Indeed during the lifetime of this journal we have witnessed the development 
of one of the most promising technological tools ever devised by man—the elec- 
tronic digital calculator. During this period the speed of digital calculators has 
increased one-hundred-thousand-fold. High-speed calculators are being used to 
translate from one language to another, to track satellites, to compose “sym- 
phonies,” to design nuclear reactors, to forecast weather, to monitor this country’s 
early warning system, to play chess, to simulate the motion of a submarine, or to 
compute fall-out patterns. This is but a small sample of the many complex tasks 
which are being performed by these devices. However, we are only in the initial 
stages of this development. It is not yet sixteen short years since the first large- 
scale automatic calculator has been placed in operation. Mathematics of Computa- 
tion will in part be devoted to the exploration of new areas of applications for high- 
speed computer devices in every field of human endeavor. 

Substantial progress has also been made in the development of numerical 
methods. For instance, until very recently it has not been feasible to obtain nu- 
merical solutions for most types of systems of ordinary and partial differential 
equations. This is no longer the case. Numerical solutions of complex systems of 
differential and integral equations—based on the methods of finite differences and 
aided by the availability of high-speed calculators—are now readily attainable. 
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2 ANNOUNCEMENT 


This constitutes a major advance in applied mathematics. The mathematical equa- 
tions which govern the solution of problems in many areas of engineering and 
theoretical physics are systems of differential equations. Supersonic and subsonic 
aerodynamics, nuclear reactor design theory, heat transfer, propagation of electro- 
magnetic and acoustic waves are but a few illustrative fields which fall in this cate- 
gory. Mathematics of Computation will be devoted to advances in numerical 
analysis with the view toward further expansion of our capabilities to obtain nu- 
merical solutions to systems of mathematical relations of any type. 

Many new mathematical disciplines are in the process of development which are 
required to advance the theory and foster the application of numerical methods. 
Information theory, the logic of automata, operational analysis, the theory of 
games, and linear programming are a few examples of such fields of current in- 
terest. These and many other mathematical developments will be required to com- 
prehend clearly and utilize fully the vast potential of modern computational de- 
vices. Mathematics of Computation will be available for the publication of advances 
in modern mathematical theories related to computation. 

In view of the phenomenally rapid progress in the development of computational 
devices and the equally impressive advances in numerical analysis, we look for- 
ward with keen interest and considerable expectation toward the future, to witness 
the exciting technological progress that will result from this major advance or 
break-through in applied mathematics. With this entire field of modern mathe- 
matics, embracing 

(1) advances in numerical analysis, 

(2) the application of numerical methods and high-speed calculator devices, 

(3) the cormputation of mathematical tables, 

(4) the velopment of new mathematical disciplines related to computation, 

(5) the theory of high-speed calculating devices and other aids to computation, 
we have associated the title, 


Mathematics of Computation. 
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Quadrature Formulas Involving Derivatives 
of the Integrand 


By Preston C. Hammer and Howard H. Wicke 


1. Introduction. In this paper we demonstrate the existence of quadrature for- 
mulas of the following types. For k an odd positive integer, 


[ = (k—1)/2 I") wo ( 
Sa) de = 2 p> @anit ya, (7° (a5) — f° (—2;)) + Rnalf)- 


For k an even positive integer, 


1 (k—2)/2 f°(0) - - 
[s@ar=2 D oat oe Uf (as) + £(—25)] + Renal S). 
1 io (22+ 1)! 
Each of formulas (1) and (2) will be shown to exist so that they are exact for 
polynomials of degree at most 4m + k for odd k, and 4m + k — 1 for even k. 
It is possible to use tables published by Hammer, Marlowe, and Stroud [2] and 
extended by Fishman [1] to obtain formulas of the form 


(3) f(x) dx = f(1) + 3 a; f"(x;) + Ra(f) 


and likewise formulas using higher derivatives. These formulas are asymmetrical 
and for some uses would be more appropriate than formula (1) or (2) 

It is considered that formula (1), (2), or (3) may be useful when the integrand 
function is represented as a variable integral for which derivatives may be easier 
to compute than the integrand function itself. It is also anticipated that these 
formulas may be used in the numerical solution of differential equations. 

Kopal [3] has devised formulas using the first derivative values. His approach 
leaves difficulties in establishing existence and reality of evaluation points and 
in some cases he has computed more than one formula of the same degree. The 
method we propose has direct connection with the established theory of orthog- 
onal polynomials which gives the existence, reality, and distinctness of the evalua- 
tion points and the positiveness of the weights a; . A linear transformation to give 
integration limits —h, h in formula (1) or (2) results in multiplying each deriva- 
tive of order n by h"™" 

We have not identified the orthogonal polynomial systems with any treated 
in detail in the literature. However, tables of the a; and x; for k = 1, 2, m = 1 
(1) 10, and k = 3, 4, m = 1 (1) 9, have been computed by G. W. Struble in 
[5], where tables for the uS") of the remainder terms are also to be found. For 
purposes of computation we give a standard type recursion formula permitting 
generation of each polynomial in a sequence from its two predecessors. 

Throughout the paper we assume that the integrand function has all-order 
derivatives appearing and that these are continuous. 





Received April 16, 1959; in revised form, July 9, 1959. This work is supported in part by 
the Office of Ordnance Research, U. S. Army. 
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4 PRESTON C. HAMMER AND HOWARD H. WICKE 


2. Reduction of the Problems. It is well known and readily established that 


1 z\ n 
(4) [ ([) g(x)(dx)"" = hd (1 — x)"g(x) dz, 
0 \vo n! Jo 


where ([ \ g(x)(dx)" denotes the result of integrating g(x) successively n 
0 / 


times over (0, x). Now our formulas (1) and (2), by choice of form, will hold for 
every odd integrand function f(x) and hence we write f(x) = fo(x) + fi(x), where 


Bt en f(x) a Sw f(x) — f(—z) 


2 


——— D end f,(2) = 


The functions fy and f; may be called the even and odd components of f, respectively. 
We now observe 


1 1 (k) k p(k) 
(5) / f(z) dx = 2/ fo(x) dx and fo" (x) = fs) + (-) Ff (2) 
La 0 





2 
Moreover, 
[ fo(x) dx = {(f) fo” (x) (dx)*** + S,(f), where 
0 
| pean y8q) 
—_____. rk 
(6) a at a for k odd 


S:(f) = 5 (k—2) /2 f?)(Q) 


| & (+1)!’ 
Hence our formulas (1) and (2) may be written, using formulas (4), (5), and (6), 
and dropping 2S,(f) from each side, as follows: 


for k even. 


9 1 m 
(7) al (1 — x)*fo (x) dx = 224; fo (aj) + Rm(f)- 


Now we only need derive formulas (7) exact for even-degree polynomials replacing 
fo(x). Then, if k is odd we have with u = 2° 


LF | (1 — x)*xP(2’) dx =if (1 — Yu)*P(u) du 


~) 


= 2 >) a; Vu;P(u;), and thus, 
j=l 


1 m 
(8) af (1 — Yu)*P(u) du = > b; P(u;), where b; = 2a; +/u;- 


On the other hand, if k is even, so that f$* (zx) is replaced by an even polynomial 
P(x’), then we have with u = z, b; = 2a;, 


1 k m 
(9) =f (= Vu) P(u) du = >> b; P(u;) 
k! Jo Vu j=1 
Now formulas (8) and (9) are in standard form to apply the theory of or- 


thogonal polynomials. The weight functions are positive-valued and appropriate. 
Hence there exists a sequence {P,,,.(u)} of orthogonal polynomials for each k = 1, 
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QUADRATURE FORMULAS INVOLVING INTEGRAND DERIVATIVES 5 


2,---, the zeros of which are the squares of the proposed evaluation points z; , 
and positive numbers 6; from which the a; are determined. The zeros are real, 
distinct, and in the interval. 

The weight functions are simple in form but we have been unable to find refer- 
ences in which this class has been treated. For k = 1 we have verified that the 
system of polynomials is not a classical one since {P%,;(u)} is not orthogonal. For 
purposes of calculating formulas, however, we state in the next section recursion 
formulas which permit sequential generation of the polynomials for each k = 1. 


3. Recursion Formula for Polynomials. The recursion formula stated here is of 
the standard sort, but its form is preferable to the one given by Szegé [4], p. 41, 
since it involves only the coefficients of P,_; and P,_, to calculate the multipliers 
B,, and C,,. This is important when explicit formulas for the multipliers have not 
been determined. We assume {P,(u)} is a system of polynomials with leading 
coefficient 1, orthogonal over the interval [a, b] with respect to the weight function 
w(u). Then we have 


(10) P,(u) = (u + B,)Pra(u) — CaPr_2(u). n 
If we define 


IIV 
bho 


b 
(11) tn = / w(u)u'P,(u) du 
we may observe u,”’ = 0,0 < k < n, and 


b b 
pn” = / w(u)u"P,(u) du = / w(u)(P,(u))* du. 
a a 


We also define c,”} as the coefficient of u’ in P,.(u). Then B, and C, are given 
by 


(n) 
‘ n—2) Mn-1 
(12) B, = Ge + =) 
Mn-1 
Pot 
‘ Y n— 
(13) C, = (n—2) * 
Mn-2 


2) 


—j) 


For our particular problem, the ui"2”, uS"7"’, and uw”): are linear combinations of 
the coefficients of P,-2 and P,_, with rational multipliers and the coefficients of 
every P,, are rational. Hence, in the absence of better formulas, starting with P» 
and P,, we can generate the remaining polynomials, in principle. In practice, of 
course, the coefficients grow very rapidly in exact rational form. 


Now Po. = 1 for every k, and, for degree one, {P;.(u)} is given by 
3 
Pe = ee — a 2. = 1,2,3,---, 
(14) Py 2 i(u) u (Qn + 3)(n + 1) n > 
and 
" 1 iti 
(15) Pion(u) =u n = 1,2,3,--- 


~ (2n + 3)(n + 1)’ 
Hence, with the recursion formula (10), successive higher-degree polynomials may 
be calculated. 
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The orthonormal sequence {Q,(«)} is determined by 


(16) P,(u) = VuS” Q,(u). 


4. Remainders for the Integration Formulas. In this section we give the re- 
mainder formulas. We will give the explicit functions to which Rolle’s Theorem 
may be applied. The method is a variation of Markoff’s for Hermite (osculating) 
interpolation (Szegé [4], p. 369.) As usual, we obtain the highest-degree poly- 
nomial H(z) for which the formula is exact and which agrees with the integrand 
function or its derivatives at all evaluation points. Then we estimate 


f(a) his H™ (x) 
and f@(xz) — H(z), and use formula (7) to find Rn x(f). (Here Ho is the even 


component of H.) 
For k odd and 2m + 1 evaluation points 0, +2, --- , +2, we require 
H(0) = f(0), HH (0) = f™(0) 

for n = 1, wee, k, H (+2;) -_ f” (+2;), H“*” (+2;) - f** (+2;),j = 1, 
- , m. Then H(z) is a polynomial of degree at most 4m + k. Let P»(2’) be the 

polynomial P,,,.(u) determining the evaluation points 2; by P,.(27) = 0, 7 = 1, 
- , m. Then let x be any number which is none of 0, +2; . Define 

f(x) — H(z) pp (2p 
- t[Pn(a z’)/? z{[F m(2 )] ° 

Then F(z) vanishes at x, 0, +2;(2m + 2 distinct points). Also, F’(z) vanishes at 

+2; and 2m + 1 other distinct points, or 4m + 1 ee Hence, applying Rolle’s 


(17) F(z) =f“(z) — H(z) - 





Theorem there exists a £ such that F“”"*”(¢,) = 0, and thus 
(4m+k+1) 
(k) («) ( ») = f ()2 [Pm (x oF 
(18) f (a) - H = a Db! . 


From equation (18) applied to z and —z, and from the continuity of f¢"***” (x) 
it follows that the even components satisfy 


, 


(b) (k) fin" (&5)2 [P..(2*)P 
(19) fo (x) — Ho (x) = (dm + 1)! : 


(4m t+ (+) 





Hence, from equation (7) by the continuity of f and the first theorem 
of the mean for integrals, we have, with u = 2’, and for some ne{—1, 1], 








4 re ) ere V/u)* P 
(20) RnaD) = "Gazal ve [Pa(u)F du, — for k odd, 


or replacing the integral by a? ; 

one (a) im 

Gm + 1)! Pe for k odd. 

For k even, we require H(x), of at most degree 4m + k — 1, such that 
H(0) = f(0), H™(0) = f™(0), 


forn = 1,---,k — 1, H“(+42;) = f®(42;), H°*?(42;) = f**(+2;), for 


(21) Ralf) = 
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for j = 1, --- , m. Then, for z not equal to one of +2;, we define 
() (k) 
(22) P(e) = f%(2) — H(2) — SE Per 
[Pn(x*)}? 


Now F’(z) vanishes at 4m distinct points, and there exists a & such that 


F“™ (£,) -— 0, 


which implies 


(4m+k) 2\32 
2: Or) — H(z) =f ev Pa 
(23) f(a) — H''(2) (4m)! ; 
The even components fo and Hy then satisfy 

(4m+k) > (2\12 
. (2) — H(z) =F) [Pale 
(24) fo (x) — Ho (x) (4m)! ; 


using the same argument used in the passage from (18) to (19). Then, as before, 
substitution in equation (7) gives, with u = 2’, and for some n¢e—1, 1] 


(4m+k) 1 pcg 
(25) Raa(f) = f ae ql ; ve [P,.(u)]} du, for k even, 
- tla/u 
or, with the integral replaced by uf? , 
f“"*® (9) A) 
(4m)! “™*? 





(26) Roa(f) = for k even. 

The number us) may be determined simply in the low-degree cases by sub- 
stituting in the numerical formula integrands x*"***" if k is odd, and x" if k is 
even. 

It may be noted that our explicit choice of H is made for convenience in each 
case. That is, formulas (1) and (2) hold with f = H and R,x. = 0. and H isa 
highest-degree polynomial for which this is true—i.e., no higher even power of x 
may be included. Moreover, in formula (6) we have 


1 1 z\ k+1 
(27) l H)(x) dx = (0) Ho (x)(dx)*" + S,(f), 
0 0 


since f”(0) = H‘ (0), for appropriate values of j. Numerical tables of a; , x, and 
the constant in the remainder have been calculated by G. W. Struble [5] for small 
values of k. 
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Tables for Use in Quadrature Formulas Involving 
Derivatives of the Integrand 


By George Struble 


1. Introduction. Tables of the a; and z; in the approximate quadrature for- 
mulas developed by Hammer and Wicke [1] 


1 (k—1)/2 (2%) 
(1) [ $2) ae = Jf (0) +> a; [f(2;) — f(—2;)], k odd 





= (21+ 1)! 
1 (k—2) /2 f° (0) | , ' 
[ f(ix)dz=2 > + Ya i (2) +f (—2;)], k even 
1 =o (22 + 1)! 
are given in table 2 for k = 1, 2, m = 1(1)10. Coefficients of related orthogonal 


polynomials are given in table 1 for k = 1, 2, m = 0(1)11. The remainder terms 
for formulas (1) and (2) are 

¥ pemeciis 
(4m + 1)! 
|s=™%9) 
{ (4m)! 


The C;.,m are listed in table 3, for the same k and m as in table 2 


| ae ; k odd 
(3) = — 


Ci.m, k even. 


2. Calculations. The polynomials P,,,.(2) were generated by orthonormaliza- 
tion of the sequence 1, x, --- , x” aa New functions 


at ,k odd. 





(4) W,(2) = | 
|( l1— aie 
k!  kla/e 
on the interval (0, 1). Their coefficients were obtained as solutions to systems of 
linear equations. The x; are the square roots of the zeros r; of Pm x(x). The num- 
bers a; are determined by 


( 


7 Am+i.k a . 
| 22; at (r;) Prnaae(1s) ; . ™ 
(5) aj =) 


,keven 





Amu, k 

\2P%, k (rj) ) P, m+1, wie (Ts) ” 

where A, is the quotient of leading coefficients of P,,, and Pm. .The values of 
Cx,m for use in the caemrueed terms were derived from 


9 2 (4m + 1)! a > a; rma : k odd 


4m +k + 2)! j= 


(4m)! = 4m 
2 6a 1s on 
2| ge 1)! 2, a | even, 


Received June 1, 1959; in revised form, August 18, 1959. The calculations and work were 
supported by the Wisconsin Alumni Research Foundation and the Office of Ordnance Research, 
U.S. Army. 


, k even 





(6) Chin - 








vere 
rch, 








Pio 


Pu 


TABLES FOR USE IN QUADRATURE FORMULAS 


© 





1.73205 08075 6888 
7.53370 80350 08842 — 2.26011 24105 0265 
30.89988 23741 8622? — 23.66207 20883 4082 + 2.68435 27159 0420 
125.09567 55614 902° — 156.54342 49462 1627 + 49.95110 87256 310z — 3.04662 58269 2541 
503.92224 31999 51z* — 878.56845 42019 08z* + 480.77283 44538 392? — 87.74164 92632 
052z + 3.36803 13248 0860 
2025.03058 73957 92° — 4532.89143 04158 Ort + 3554.73528 18164 Ox? — 1152.71421 91499 
2x? + 138.22772 15120 70z — 3.65996 50721 1686 
8126.60217 35502 z* — 22225.89230 73142° + 22794.37521 94452* — 10831.68147 28542? + 
2371.48444 08556 z? — 202.49633 19667 42 + 3.92935 16425 271 
32585.34917 36427 — 1 05329.15543 67z* + 1 33708.35590 972° — 84369.53936 860zc* + 
27562.57790 4932* — 4390.47486 255827 + 281.55199 516262 — 4.18076 36516 86 
1 30585.72574 2* — 4 87145.05832 x? + 7 37945.40041 2° — 582772.09740 r> + 2 56072.42057 
z* — 61782.59435 92° + 7521.30300 462? — 376.33313 306z + 4.41740 68709 
= § 23121.3992° — 22 12241.352° + 38 96351.9527 — 37 02822.37r* + 20 55853.62z5 — 

6 74493.5592* + 1 25947.1282* — 12137.51852? + 487.72367 9x — 4.63163 193 
= 20 95022.52'° — 99 04599.42° + 198 92582.2° — 221 24981.77 + 148 91824.2° — 

62 25087.32° + 15 95661.1z* — 2 38415.82z* + 18678.058z* — 616.56109z + 4.85521 90 
= 83 88194.2"! — 438 42170.r'° + 989 11790.2° — 1261 15400.c* + 999 31000.27 — 
509 73540.2° + 167 70010.72° — 34 68084.2* + 4 25177.7z* — 27649.292? + 763.60752 — 
5.05931 2 


ll i] ll 


1.73205 08075 6888 

12.70977 81860 4492 — 1.27097 78186 0449 

59.88757 23920 5992? — 29.17599 68063 882z + 1.20652 61837 2282 

259.32961 55268 10z* — 242.98324 64964 47z* + 50.37043 53148 982r — 

1.18192 23907 2002 

1087.36971 21542 1lz* — 1531.67966 14027 7 x* + 633.73798 31134 97z? — 

76.31728 62592 1962 + 1.16909 13241 5856 

4487 .70266 55197 82° — 8483.52105 08986 2z* + 5405.81028 77031 Or? — 

1335.91407 05265 8x? + 106.99341 52181 20z — 1.16124 14205 1600 

= 18355.66707 84232° — 43644.30820 73452° + 37942.62659 31132* — 14721.14915 188423 + 
2472.29781 31355 x? — 142.37405 03369 32 + 1.15594 79227 687 

74659.95872 21927 — 2 14151.18872 092° + 2 36965.50353 3125 — 1 27161.63787 63r* + 

34058.68565 0872? — 4184.54516 34292? + 182.43804 37977z — 1.15213 67777 11 

= 3 02544.85398 62% — 10 16908.12494 z7 + 13 69647.20094 z*° — 9 45055.96574 Or*® + 
3 54167.19148 2z* — 70308.79999 50 x* + 6633.07482 9847? — 227.16794 49007 + 

1.14926 09079 5 

12 22819.38562° — 47 14584.9219z* + 74 930684.8783z7 — 63 42588.4460z° + 

30 83282.81842° — 8 64946.912662* + 1 33312.130542* — 9996.98758 liz? + 

276.54933 8282 — 1.14701 28163 

= 49 33033.9142'° — 214 63128.26r° + 393 48240.15 z* — 395 46800.17z7 + 237 46979.27z° 
— 87 07972.887z> + 19 12893.1892* — 2 36576.7122z* + 14474.00552 xz? — 330.57026 87x 

+ 1.14520 6765 

198 72572.1z"! — 963 24095.72'° + 2001 94335.2° — 2332 57625.c* + 1672 30868.27 — 

762 13570.12§ + 220 33307.37° — 39 12962.352* + 3 98072.797z* — 20280.6219 x? + 

389.22470 8x — 1.14373 626 


ll 
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TABLE 2 
k=1 
m —E _ 
xj aj 
1 .5477225575 05166 . 3042903097 25092 
1 .3721456511 38755 .2998694151 77678 
2 .7920059217 60865 0695342880 47154 
1 . 2820900214 17742 2625168155 42158 
2 .6268601608 87658 1158516504 33743 
3 8825310979 63249 0226513358 72470 
4 1 . 2273288824 80793 . 2283023164 62706 
2 .5140806061 43694 .1313690234 71769 
3 -7563681074 02448 0510204966 76312 
4 9248839673 29570 0093443207 77566 
5 “- .1904841804 26525 . 2005148719 02275 
2 | 4344049427 27343 -1330163257 25257 
3 | 6535558761 15115 0693217404 15712 
4 | . 8292980447 85322 0254519669 74067 
5 .9479286287 01912 0045107921 65497 
i 2 | . 1639811690 9936 .1782134815 8969 
| 92 | . 3756528865 1806 . 1290409387 9439 
_ = .5724366400 9800 .0790426562 2466 
whee Y .7417866536 8859 0393058004 7950 
| 5 | 8740336979 4142 | 0139874045 8116 
| 6 9618129567 5531 0024324030 3689 
| | 
1.2 | . 1439912886 45 . 1601387035 18 
am .3307098113 32 1229802576 16 
} 3 | .5079816905 42 0833391819 77 
| 4 6669327097 06 .0490821180 75 
| 5 | . 8006064071 09 0237352491 28 
| 6 | .9033361608 21 | .0082842932 84 
oe I .9708112198 58 - 0014235185 93 
| | | 
bi 4 1283703842 . 1452769503 
2 . 2952852684 1163542537 
. | .4559661731 0844453524 
el 6038699343 .0553533943 
5 .7338622714 .0318100298 
6 | 8415994425 0151045544 
oy 9235314840 0052043396 
8 | 9769705854 0008871299 
1 115823631 . 132878407 
2 . 266673730 -109816541 
3 413291816 083725810 
4 550684038 059056710 
5 674964515 037887335 
6 .782771237 021412447 
7 871235919 .010039513 
- 938022434 003428033 
a 981369607 .000581047 
} | | 
i | 1 | . 10552258 12239571 
2 | 24309589 10363706 
ma .37773918 | .08198747 
| 4 50552941 06098113 
| 5 62345737 04220687 
| 6 .72888457 02667176 
a . 81949694 .01489473 
= 89332330 .00691990 
ge 94876470 00234732 
| 10 


- 98461980 


-00039623 
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TABLE 2—Continued 





k=2 


xj | aj 








~I 


10 


nue 


wn 





wm Ode 


or Whe 


or wndoe 


a 


NOOO Whe 


ore COON 


— 
Sc 


on 


CBN Wh 





SCUVanNouarwnre 


_ 


-3162277660 16838 


- 2136036212 56443 
6644945298 23701 


- 1638297723 39671 
-5118494949 80572 
- 8050693913 12212 


. 1336140770 13473 
-4169305407 03968 
-§739890075 18818 
- 8733068456 92833 


- 1131072568 64736 
- 3520392516 58368 
. 5769763210 36054 
- 7683849647 09291 
-9112358854 64930 


-0982014632 7176 
- 3048169927 4399 
- 5034255107 0724 
-6810857911 4211 
-8274681084 1897 
-9344124389 0913 


-0868450485 199 
- 2688827298 199 
-4461187708 020 


- 6095094010 
- 7517394338 839 
86669856 


75 800 


“9495850305 493 


-0778891768 89 
- 2406041327 18 
4003534653 


-5505429051 14 
- 6858231067 44 
8015939596 43 
- 8940032842 24 
- 9600499057 83 


-0706369010 
- 2177597386 
- 3630238378 
- 5014536828 
- 6290237349 
- 7422126043 
-8379675527 
-9137418873 
-9675684295 


ooh 


CwrOows! 


-0646394150 
- 1989137370 
3320234405 


-4601078372 
-5800614933 
-6891397684 
- 7849102517 
-8652666620 


- 9284600004 
- 9731500562 





. 1666666666 66667 
. 1437779500 85321 
-0228887165 81346 


- 1224797344 91539 
0395183540 10279 
.0046685781 64849 


- 1060360979 05865 
.0469644870 94641 
-0123633215 15102 
-0013027601 51060 


-0933196598 68176 
0494739624 16079 
-0188805340 81087 
-0045407258 01483 
-0004517844 99852 


-0832757144 1884 
-0495578221 5891 
-0234294328 9172 
-0083180426 8343 
-0019025149 7640 
-0001831395 3766 


-0751675980 885 
-0484760818 692 


958 602 


-0263277' 

-0117357442 186 
0039910960 716 
0008849203 899 
-0000834301 748 


0684940996 00 
0468585232 40 
-0280343978 67 
-0145005204 06 
0062300015 19 
-0020600871 63 
-0004474372 29 
-0000415997 26 


0629090570 4 
-0450255142 4 
.0289225611 3 
.0166061935 8 
-0083212591 4 
-0034868683 1 
-0011308362 0 
-0002420894 9 
-0000222879 4 


0581677528 
-0431420181 
-0292595373 
-0181452843 
-0101353888 
0049671598 
-0020445253 
-0006538257 
-0001385086 
-0000126607 
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TABLE 3* 
b= 1 k=2 

ae ee ie Meee ae Gin ; 
1 1.7619 xX 10° 1 6.1905 xX 10-3 
2 1.0473 K 10° 2 2.7882 xX 10> 
3 6.3902 « 10-5 3 1.4869 x 10-5 
4 3.9380 x 10-6 4 8.4576 XK 107 
5 2.4386 x 1077 5 4.9654 x 10° 
6 1.5142 x 10° 6 2.9680 x 10-° 
9.4185 x 107° 7 1.7940 x 107° 
8 6.3112 <x i 8 1.0922 x 10-™ 
9 we x io 9 6.8464 x 10-"% 
10 —3.3 <x 190° 10 —2.8 x 


10-3 








* These figures include round-off errors. 


and thus show effects of round-off errors in a; and 2; as well as the theoretical 
remainder. Since the matrix of the linear system solved to find the coefficients of 
Pn.«(x) is ill-conditioned, the accuracy of the figures decreases rapidly with in- 
creasing m. Figures are kept to the upper limit of accuracy. If results are rounded 
to one or two fewer significant figures than are carried in the table of x; , there is 
no doubt of the accuracy of the digits kept. 

The formulas were derived in [1]. The calculations were carried out in the 
Numerical Analysis Laboratory of the University of Wisconsin on the IBM 650. 
An 18-digit floating-decimal interpretive routine by Eugene Albright and a linear 
systems solver by Gerald Thorne were used. I am indebted to William Kammerer 
for several helpful suggestions. - 


University of Wisconsin 
Madison, Wisconsin 


1. P.C. Hammer & H. H. Wicks, “Quadrature formulas involving derivatives of the inte- 
grand,’’ Math. Comp. (MTAC), v. 14, 1960, p. 3-7. 
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Numerical Quadrature Over a Rectangular 
Domain in Two or More Dimensions 


Part 1. Quadrature over a Square, Using up to Sixteen Equally Spaced Points 
By J. C. P. Miller 


1. Introduction. Except for a section of a paper by Bickley [1] which gives some 
of the results that follow, there seems to be very little in print concerning numerical 
quadrature over rectangular domains in two or more dimensions. This is perhaps 
because numerical evaluations can be readily made by using one-dimensional 
formulas on each variable in turn—such formulas as the Euler-Maclaurin formula, 
Gregory’s formula, Simpson’s rule, Newton’s three-eighths rule, or Gauss’s formu- 
las, to name a few. This restriction to “products” of one-dimensional formulas 
limits the field unduly, and may lead to an excessively large amount of work. 
Thus, in n dimensions, use of Gauss’s 3-point formula involves 3” points, whereas 
comparable accuracy may be obtained with about 2n? points. 

In this first note we explore some of the simple possibilities corresponding to 
Simpson’s rule and the three-eighths rule, applied to 9 or 16 points equally spaced 
over a square, the corner points being included. 


2. Integration over a Square. We restrict the domain of integration to be a 
square; this covers any rectangular domain by change of scale. We take the center 
of the square as origin of coordinates, and in the first place take the side of the 
square to be 2h; and the integral as 


h h 
(2.1) 1=[ [f(a u) dz ay. 
h h 


We assume that f(z, y) can be expanded as far as we need in a Taylor series in 
x and y. In other words, we suppose that f(z, y) can be represented adequately by 
a polynomial in z and y, with an error term which we shall suppose may be esti- 
mated by considering a few of the more significant neglected terms. We shall not 
give an accurate error analysis. 

It is clear that polynomial terms involving an odd power of either variable 
will contribute nothing to the integral; we shall also group ordinates in sets such 
that the total contribution to their sum is zero for such terms involving an odd 
power of either variable. For example, the points (h, 0), (—h, 0), (0, kh), (0, —h) 
form one such group. With this grouping we can then eliminate from considera- 
tion all terms of the Taylor expansion which do not involve an even power of both 
variables. 


3. Method of Derivation. We seek approximate formulas of the form 
h h 
(3.1) I = | [ f(x, y) dx dy = >> A,,.(rh, sh) , 
-h h 


Received September 18, 1959. 
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in which A,,, = Az, = Ajrj,,5, and which will be exact for appropriate poly- 
nomial functions f(z, y) of sufficiently low degree. 

Two approaches, not entirely distinct, may be used: 

(i) We may give f(z, y) special forms and choose coefficients A,,, to fit these 
exactly. Such forms are 1, z* (or the equivalent symmetrical form 2? + 7), z‘*, 
xy’, ete 

(ii) We may expand f(z, y) as a Taylor series and evaluate J both by means 
of the integral and by means of the sum, and equate coefficients. 

We shall give an example of each approach, in order to bring out an interesting 
point concerning the proper choice of special forms. 


4. Nine-point Formulas (i). Consider the nine points, grouped as indicated 
by the semi-colons: (0, 0); (h, 0), (—h, 0), (0, h), (0, —h); (h, h), (—h, h), 
(h, —h), (—h, —h); and the four special functions f(z, y) = 1, 2’, x‘, 2°y’ with 
values at the nine points as follows 


1 Vee ey ee eee 


1 1 
(4.1) 1 Frit | o.3.3,- 1.5.8.8 0 
1 1 








1 101] | 


i -_- oOo — 


The values of 1/4h’ are 
(4.2) 1, 


Hence, we wish to have 


ool 
ode 
a 


(Aoo = 4Aio + 4A = 1 

‘ 2Aio + 4Ai, _ $ 
(4.3) 2A1,.0 + 4A = $ 
4Aj4 _ $ 


We note at once that two of the equations are incompatible. Since we have 4 
equations for 3 unknowns, we can still solve them if we discard the third (since 
the correct result for x” must have precedence over the correct result for x). We 
obtain 

Ao,o _ $ Ao _ $ Ain - 3 


which is precisely the Simpson’s rule “‘product,” see Bickley [1], eq. (22), with 
multipliers 


es eee a ae 

it 4 1/1428 
(A) [4 16 4 | 
1 41] 


to give [/4h’. The main error term clearly concerns the coefficients of the terms 


in x‘ and y/‘ and is 
h’ (af + $3) 
180 ax! oy 0,0. 








oly- 
1ese 
2, 
ans 
ting 
ated 


h), 


with 


ive 4 
since 
. We 


with 
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5. Nine-point Formulas (ii). We now consider the same problem by expanding 
and equating coefficients. We have, using f,., for f(rh, sh), and fo simply for f(0, 0), 


2 a2 » 22 4 
- x 9 fo y ofo x* d'fo 6ry Ifo ¥ ahe 
x P = — ee Se «.e eee 
(5.1) f(z,y) =fo Sia + Daye + Tart 4! axtaye * May + 





whence 
— Ifo , Ifo dfs , 10 df . d'fe ‘ 
— [s +3 Ge id jaa 3 ni (ss + 3 aay t ch) de | 


Following Bickley [1], we write this concisely in terms of 





ee a 
5. Yal+— onl 
(5.3) aa tae D ~ axtay?” 


Then, adding further terms for later use, 


2 2 6 
1/4h? = fo + + Vf + = (V*fo + $D'fo) + (V°fo + 4D°V'fo) 
(5.4) = . 8 ‘ 
aa = (V*fo + 8D'V'fo + 48D%fo) + --- 
We have likewise 


fo.o = fo 
_ 2 4 
fue + for + t10 + faa =Afo + 5 RV + = (Wife — 20) 


+ 2 fe — 39'V'fo) + 2 Te — 49'V'f. + 2D%fo) + --- 


9. 9 4 + 
fia +t faa tha t+ faa = 4fot+ 5 ivf + ry h'(V'fo + 4D fo) 





+4 = WV fo + 12D'Vfo) + = W(Vfo + 24D'Vfo + 160%) + --- 


From these relations we obtain, for terms to A* in 1/4h° 


Aoo + 4Aj9 + 4A31 = 1 
2Ai0+ 4A, = \ 
2Ai0 + 4A _ 5 

| = 4Aio a 16Aj, - i's 


(5.6) 


The first three equations are as in (4.3), and the second and third remain in- 
compatible. Neglecting the third we can solve to give (see Bickley [1], eq. (20)) 
the multipliers for 1/4h’; 


(B) 16 88 16 The main error term is 180 V'fo. 
1 4 
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The difference in the fourth equation between (4.3) and (5.6) leads to different 
formulas. The cause is reflected in the error term. In the second case we have 
ignored a term in V‘fy and removed the remaining term in fo ; in 4 we have ignored 


a term in 
O'fo 4 =*) 
Ox oy 


bees . 4 . - 
and removed the remaining term in D fo. The error term in the product-Simpson 
formula can be expressed alternatively in the form 


4 
7 (V*fo — 2D°fo). 

It is perhaps a moot point which formula will give better results, unless it is 
known that f(z, y) is exactly, or almost a harmonic function. It is perhaps worth 
noting that for the method of 4 to yield the formula (B) it is necessary only to 
take a harmonic function instead of 2’y’ in deriving the last equation. For ex- 
ample, f(x, y) = «* — 62°y’ + y' gives values 


p—a 3 ei 

(5.7) | 10 1| and J/4h? = — 4h‘ 
|'-4 1 —4 

and the fourth equation 

(5.8) —4Aio + 16A,4 = is. 


6. Five-point and Other Formulas. The impossibility of removing all of the 
h* terms in the error in J/4h’ by using 9 points, suggests the possibility of using 
fewer points, in fact five or eight, while still retaining an error of order h*. There 
are three possibilities, using the first two equations only of (4.3) or (5.6). 

(i) Take Ao.» = 0; this yields multipliers 


-———— 








| 
|}-1 4 -1|/+12 
(C) | 40 4 | with main error y}gh*(V*fo — 22D*fo). 
—-1 4 —-1 
(ii) Take A; = 0, giving 
i. =< coe 
01 0 | + 6 
(D) ee a with main error ¢ish*(V*fo — 7D*fo) 
eo i 6 


14 
(E) | 0 





f 


nt 
ve 


ed 


on 


; is 
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to 
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“fo a, 
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The first of these has little to commend it, with its large D*f, error term, and 
its negative multipliers. The other two (see Bickley [1], eqs. (23, 24)) are reason- 
ably good. Over large areas, with many contiguous squares of area 4h’, they aver- 
age respectively 3 and 2 points per square. With a square of side 6h, for example, 
they give 

















010101 O|+6 11020201\+2 
1223331 10 80808 0 
es:e80e023 6 1204040 2 
(F) 11222221 and |0 808 08 0 
020202 0 1204040 2] 
12232232321 080808 0| 
0101010 }1 02020 1) 











The first of these has interior points with multipliers 2 (or zero) only. We note 
also that combination in proportions 34:75 gives precisely the formula (B). One 
further formula is obtained by combining (D) and (E) in equal proportions: 


1 + 24 


(G) 2 ; with main error zigh*(V‘fo + 4D°fo) 
Hi) 
and with simple multipliers and a good error term. 


7. 16-point Formulas. These follow exactly the pattern for nine-point formulas; 
we merely quote results. 


The method of 4, with the same special functions, and with similar neglect of 
the incompatible third equation gives for I/9h°, where 


3h/2 3h/2 
= [ [ f(a, y) da, dy, 
3h/2 3h/2 


the multipliers 


& 4 + 64 


Te = 
lig 9 9 : 

(A’) : ; : . with main error term ¢gh'(V‘fo — 2D'fo). 
13 3 1 





The formula corresponding to (B) of 5 is 





7 13 13 7 | + 320 


‘ 13 47 47 13 “ - ~e 
(B’) 13 47 47 13 with main error term ;yh'V fo. 
i? 2 ae? 





8. Twelve-point Formulas and Others. As in 6 we can use the first two equa- 
tions after setting one of the coefficients A;,; , A1,3 , or A3,3 equal to zero. This yields 
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twelve-point formulas 
(i) with Aii = 0: 











—2 3 3 -2|+16 
(C’) - - : : | with main error term 5h‘(V‘fo — 47D'fo) 
-2 3 3 -2| 
(ii) with A33 = 0: 
01 21 0/+ 16 
, »2 2 2b] ote aoe 
(D’) 1221 with main error term gh (Vfo — 7D fo) 
011 0) 


‘ 
(E’) lo 
F 





Also, as in 6, we may combine the last two in the ratio 4:4, to give 


| 


122 1\+48 
f, i277 2) 103g = 
(G’) 277 2 with main error term @yh'(Vifo — 4D fo). 
a a aa | | r 





Finally, we write out multipliers over a square of side 6h for the 12-point formu- 


las (ii) and (iii) above (which average respectively 8 and 5 points per square of 
side 3h) 














0110411 0);+16 1002 0 0 1|+16 
ha @2 3 2 i 03 3 0 3 3 0 
iS) S23 Bd 033 0 3 3 0 

(F’) o2202 2 0 1004001 
fee aead i 03 3 03 8 0 
oS 2 = a a 03 3 03 3 0 
011011 0 1002001 

















Note once again the multipliers 2 or 0 only in the interior in the first case. 


9. Numerical Illustrations and Comments. We consider the application of the 


formulas to two examples: (i) f(z, y) = cos x cos y;and (ii) f(z, y) = sin z sinh y, 
a harmonic function. 


1 1 
(r= | [e082 008 y dz dy = 4 sin* 1 
1 1 


whence 17 = 0.708073. 








NUMERICAL QUADRATURE IN TWO OR MORE DIMENSIONS 19 


The series (5.4) gives 


fo 1.000000 

V'fo/3! — .333333 
V'fo/5! + 33333 

4 D'fo/3-5! + 11111 
V‘fo/7! _ 1587 

4 D'V'f./7! _ 1587 
V*fo/9! + 44 

8 D'V*f./9! + 88 
16 D'fo/5-9! + 9 
Sum 0.708078 


This indicates the need for another term, and the alternation in sign, in this case, 
of terms with successive powers of h’. 

Table I shows the results of applying various formulas, with various values of 
h. Results are given to 5 decimals (based on calculations with 6-decimal function 
values). The column «¢ gives the error (formula—true value) and column C gives 
the computed value of the leading correction term, as listed in earlier paragraphs; 
« and C are in units of the 5th decimal. 














TaBLeE I 
h=1 h =1/2 h = 1/3 
Formula | ; 
| Reutt | « | Cc Result | « | C | Result ‘ 
(A) 0.71701 +894 |—1111 (0.70858 | +51 | —53 (0.70817 +10 
(B) | .72641 | +1834 |—2222 | .70909 | +102 | —107 | | 





(C) | .62309 | —8498 '+10000 | .70345 | —462 | +481 | | 
(D) | 169353 | —1454 |-+1667 | -76730 | —77 | +80 0.70792 | —15 
(B) | .76308 | +5501 |—6667 | .71114 | +307 | —320 

| 




















(D’) | .70175 | —632 |+741 .70773 | —34 | +36 


(G) | .72876 | +2069 |-2500 | .70922 | +115 | —120 
h = 2/3 ares hk = 1/3 e ; 
Formula 
Result € Cc Result € Cc 
(A’) |0.71199 | +392 \_494 0.70830 | +23 | —24 | 
(B’) | .71608 | +801 |-988 | .70852| +45 | —48 
(c’) | .61988 —8819 |+10617 | .70319 | —488 | +511 


























(E’) -74269 | +3462 |—4198 | .71000 | +193 — 202 
(G’) - 71540 +733 |—905 | .70849 | +42 —44 
Taste II 
h = 0.6 h=04 

Formula —| Formula a 

Result € | c | Result € Cc 
(A) | 0.129414 | +187 | 186 | (A’) 0.129310 | +83 | —83 
(B) 129227 | | | (BY 1292297 | 0 | 
(D) . 129879 | +083 | —652 | (D’) .129517 | +290 | —290 
(E) | a | —745 | +745 | (B’) | .128689| —538 | +538 
(G) | :120181 | —46 | “447 | (G) | [199081 | 414 | —14 
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We note that the correction estimate C is numerically larger than the actual 
error ¢ in each case; this is due to the sign alternation mentioned earlier. We see 
also that, even with h = 1, C is quite a reasonable estimate of the size of «. 

Formulas (C), (C’), as expected, are not very good. Formulas (E), (E’), are 
also poor, since V‘fy and D*f, have the same sign. The “‘product-Simpson”’ rule 
gives the best results—best of all with h = 4. Forraulas (D) and (D’), as exhibited 
in (F), (F’) for a square of side 6h, suggest that (D) gives better results than 
(D’) when the same points are used. 

In practice, it is useful to use two formulas, and to compare results to give an 
estimate of the possible error. For this purpose (A) and (D) or (A’) and (D’) 
seem suitable pairs, unless it is known that f is harmonic or V‘fy is expected to be 
small compared with D*fy . 


1,2 1,2 
(ii) I =-/ i sin x sinh y dx dy = (1 — cos 1.2)(cosh 1.2 — 1) = 0.516908. 
0 0 


Approximations to 7/4 = 0.1292271 are listed; all 6 working decimals are shown 
in Table II. In this case, for a harmonic integrand, the superiority of (B) and 
(B’) is evident; (G’) is also good. In (B) and (B’), since V’fo = 0, the error terms 
are of order h’ instead of h*. The quadrature of a harmonic function will be con- 
sidered further in a later note. 

I am glad to acknowledge help received with numerical calculations from Dr. J. 
W. Weench, Jr. in Washington, D. C., and from W. R. Rosenkrantz at the Univer- 
sity of Illinois; I am also grateful for facilities placed at my disposal at the Digital 
Computer Laboratory of the University of Illinois. 


The University Mathematical Laboratory 
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The Digital Computer Laboratory 

University of Illinois, Urbana, Illinois. . 


1. W. G. Bickuey, ‘Finite difference formulae for the square lattice,’ Quart. Jn. Mech. 
and Appl. Math., v. 1, 1948, p. 35-42. 
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Numerical Integration Formulas of Degree Two 


By A. H. Stroud 


1. Introduction. Here we discuss numerical integration formulas of the form 


[ H2)w(e) dz = > a;f(v,) 


where F# is a region in n-dimensional, real, euclidean space; x = (2, x2, -°-- , Zn); 
the a; are constants; and the »; are points in the space. Most previous authors 
have given formulas for special regions (for a bibliography see [4]). Thacher [7] 
has given a method for constructing formulas of degree 2 with n + 1 points for 
general regions and of degree 3 with 2n points for certain symmetric regions; with 
his method, however, each region must also be treated separately. Our main re- 
sults are to obtain specific formulas of degree 2 with n + 1 points for a general 
region satisfying a certain condition of non-degeneracy, and to show that for 
these regions such formulas cannot be obtained with fewer points. We also give a 
specific 2n point formula of degree 3 for a general centrally symmetric region. 
These results are a generalization of those of Georgiev [1, 2, 3] who has obtained 
similar results (but gives no general formulas) for n = 2, 3 with w(x) = 1. Our 
results are obtained by a different method which was developed without knowledge 
of Georgiev’s work. 


2. Formulas of degree 2. We assume at first that an integration formula of 
degree 2 for R with respect to w(x) can be obtained with n + 1 points 


Vs = (Pn, °° 5 Pin)s ¢=0,1,- 


» n, 
Then the equations 
ao + a +--+ +4, = &% 
(1) om Fang tess + Gnrng = Cj 
AovosVor + AvijVie + +++ 1 OnMajVnk = Cie j,& = 1,2,--- a 
must be solved for both the a; and the v; , where 
Co = i w(x) dz, Coj = / x; w(x) dz, Cy = / 2; x, w(x) dz. 
R R R 
We begin by writing (1) as the matrix equation 
(2) U"AU =C 
where 
} Vor *** Von am 0 --- O Co Co °°* Con 
a lm -**) Mn aes ye. 6 @ 7 
, Yat *°* Pan ». &. +s ey — a =" Ge 


and where we assume 0 < co < «© andO <|detC| < <. 


Received March 18, 1958; in revised form August 4, 1959. This work was supported in 
part by the Office of Ordnance Research, U. 8. Army and in part by the Wisconsin Alumni 
Research Foundation. 
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Since C is non-singular we can find a matrix T such that 
(3) T°U°AUT = T°CT = @E 


where £ is a diagonal matrix with elements +1. The method for finding T is well 
known (see [5], p. 56); we illustrate it using n = 3. 
Since co ¥ 0 we define th; = —co;/co, i = 1, 2, 3, and form 


, Co 0 0 0 
to too tog 10 
r i. 1 0 0 6 an eter wl” ci?” cf) chy 
he Sosa Samm of cP of | 
so 2 4 

O cis cy cy 
Now if ci)” = 0 some cj; ¥ 0 since det C ¥ 0. Assuming ci}? ~ 0 we form 
1000 Co 0 0 0 

010 4 es O heir’ + h'edy fy + hey) cf + heiy 

ian 0h10 C= 12 Gt: = (1) (1) (1) (1) 
. 00 4 0 Cie + her C22 C23 
O ety + hess c33 cy 


and choose h so that cf?’ = 2he{?? + hci? ¥ 0; if cl?” ¥ 0 we take h = Oso that 
ci) = c{1’. In this way we are assured that the element in the 1, 1 position is ~0. 


Similarly we may find matrices 7; , 7, and 7’; such that 





Co 0 0 0 | 
0 0 0 Pa ss! 
Cy = Te TC: T3 Ts = » | *@ab=9 9 & 9 | 
0 0 ci C33 0 = | 
(2) (2) 0 0 0 esa) J 
0 0 C22 C33 | 


where c¢3;’ and c3; are <0. Defining 7’. as the diagonal matrix 
[1, [co/| eh? |}? , [co/| eS |}? , feo/| eS? {1') 
we have finally T = 7,7T2T3T,TsT-¢ . 

We can assume F has the form [1, 1, --- , 1, —1, --- , —1] since any other 
arrangement of +1’s and —1’s can be put into this form by a suitable interchange 
of the rows of UT and the corresponding columns of T* U". If C is positive definite 
(for example if w(x) is of constant sign on R) EF will be the identity. It should 


be noted that the first element of F will always be positive. 
In the following we write 


1 gu s+ Eon 1 vo cw 4 4 TO01 es ~ 
UT = l fu Ka ° Ein ns 1 vu os - 0 Til . <2 P Tin ; 
1 £n1 eee Enn 1 Val eee Nes ® Trl ie naa 


Because UT is non-singular and FE = E we easily obtain from (3) 


(UT)E(UT)* = A. 





Ir 


(4 


Ww 


of 


ell 


hat 
~(), 


L_-___ S 


her 
nge 
nite 
yuld 





NUMERICAL INTEGRATION FORMULAS OF DEGREE 2 23 


In terms of the £; this equation is 


(4) Lt Enka to + bin in — Bi Batt — ++ — bin bn = By 


i,j, =0,1,---,n. 
where p + 1,0 S p &S ni, is the number of +1’s in E. We discuss the solution 
of (4); the »; are obtained from the &; by »;; = 70; + Eatiy + --- + Emtay, i = 
0,1,---,n,j = 1, ---,n, where 

1 ro ++ Ton | 
yt Oru -': Tin 


, ’ 
0 Tnl ae See 





We are only interested in real solutions of (1) and therefore precisely n — p + 1 
of the a; must be negative by Sylvester’s “law of inertia” ((5], p. 56). If E is the 
identity (p = n) clearly we must have 0 < a; < «@ ; if p < n the only condition 
for the a; is that they be non-zero. 

Table 1 gives a particular solution of (4); we have assumed ap , --- , @,_» nega- 
tive and a,p41,°°*, @, positive. In the places where a double sign occurs we 
mean to use the lower sign for the last n — p components of each vector and the 
upper sign for the first p components. Each £; is real. 


TABLE 1 











sie aibs. it tet 1/2 1/2 12 
akc! Co(Co — Go — a; — G2) + 00a = ay ) 
+(co — ay) — ai)a2 (co — do)(co — ao — a) Co — 





a 1/2 1/2 
&=(0,0,--0 =- 2-2 * ao ) 
sco — a)ai Co — Ao 


0(co — ap — --* — Gn) |? 
0, = 
(co -— ayo o+> = Gn-3)An-2 
“ +a, 1/2 a +c eo 1/2 e tay a 
> “1 (co — a — ai)(co — ao — a1 — a2) | "| (co — ae)(co— ao — a) | * | co — 
s00(co — do — --* — Gn) i? +00On_2 ug 
gn-1 = 5 — = °°? 
(Co — Ap — *** — Gn-2)Gn-1 (Co — Go — -** — Gn_s)(Co — Gp — -** — Gn-2) 
; (co — ao — a1)(co — Go — a — Ge) f (co — ao)(co — ao — a) : Co — A 
+004n-1 1/2 +09 Gn» 1/2 
En = + ’ + ates 
(co — Go — +++ — Gn-2)Qn (co — Go — *** — Gn-z)(Co — Go — *** — Gn-2) 


+0942 1/2 +094; 1/2 +a 1/2 
+o, — e ‘ +- a an ae 
(co — @o — a;)(co — ao — a; — 2) (co — ao)(co — ao — a) Co — a 
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From a particular solution £;; of (4) other solutions may be obtained as follows. 


If 
} 0 --- O 
ga }O ou vs om 
- ta. .*2* dee 


is a cogredient automorph of E, that is if SES* = E, then 

bij = Eaoay + +++ + Een; 
is also a solution. If Q is an arbitrary skew matrix of order n + 1, with first row 
and column entirely zero, such that det (EF + Q)(E — Q) # 0, then 

S = (E+ Q)'(E — Q) 


is a cogredient automorph of FE (see [5], p. 65) of the above form. If E is the identity 
S is orthogonal. We remark that in this latter case (4) determines the distances 
d(é;, 0) and d(é, ’ £5), ‘iJ ” 0, 1, <> 5 i # ds 


dé; ’ 0) sag [(co a a,)/a,} dé; ’ £;) = [co(a; + a;)/aa;}' 
The formulas discussed above are minimal; that is, similar formulas cannot 


be obtained with fewer points. For if a formula could be obtained with m + 1 
points »;,7 = 0, 1,---,m, m < n, then equation (2) would still hold, where C 


is the same as before and 
F M1 2° ee vs S i sania ° | 
eT = t ms +++ ian Sf theca e @ ' 
tare Ta gee 


that is, U is a rectangular matrix. Since U and A have rank at most m + 1, then 
U* AU has rank at most m + 1 and therefore det (U"AU) = 0. By assumption 
det C ¥ 0 and thus (2) cannot hold for m < n. 


3. Formulas of degree 3 for centrally symmetric regions. We assume RF to be 
centrally symmetric with respect to the origin; then if x is in R, —z is also in R. 
Let us further assume w(—2x) = w(x) for x in R. Then 


/ z;w(x) dx = | 2,x;7, w(x) dx = 0, 6,j,R =m 1,-+> 
R R 


We may obtain an integration formula of degree 3 for R with respect to w(x) with 
2n points as follows. Take the points to be »;, —v;,7 = 1, --- , n, and take », 
—»v, to have common weight a,. Any 2n points chosen in this way integrate 
exactly the monomials 2; , x, ja, with respect to w(x) over R. In addition we must 
solve 


a + a soe my = 3 


Ayr ;Vix + AnvejVee + ++ + Anvnj¥nk = 3Cje j,k = 1, +++, 0. 
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The second of these may be written as the matrix equation (2) where now 


— —— Vin aq 0 eee 0 G1 Ce -*- Cin 

U — oo -** fe 4 Om --- 0 C Ll] Giz Con ++ Com 
j = Lt = =>=_— 

* Vn2 ne? Van . 0 -oe a, Cin Con eee Con 


and where we assume — © < @ < ~ and0 < |detC| < o. 


We solve this equation by a method similar to that of the preceding section. 
We find a non-singular matrix T such that 


T°U"AUT =T'CT=E 
where E is diagonal with elements +1. Again it is convenient to assume 
E = [l,---,1, —1,---, —1] 


where the first p elements are +1,0 < p S n. Now writing 


bn £12 o" Ein M1 «=~Vi2 eins * 1 fen Ti2 *** Tin 
UT = fo - Eon |_| va - Von T1 Te *** Tan 
Ent Ene Mg ding Ean Vni —‘Vn2 ala bi al B Tn2 pie: Tnn 
the £;; may be solved for in terms of the a; . This gives 
- Rw 7 
(5) takat --- + Sip bin — Exper Eps — °** — Sin bin = ae tae I,--+ 
precisely n — p of the a; must be negative in order that the £; be real. 
If a, --- , a, are positive and a,4;, --- , a, negative a particular solution of 
(5) is 
& = (0,---,0, V1/|a;|,0, --- , 0) i=1,---,n 
where the 7th component of £; is non-zero. If S = (¢;;) is any cogredient auto- 


morph of FE then “a = £101; + --- + mon; is also a solution of (5). If EZ is the 
identity, that is, C is positive-definite, the solutions of (5) correspond to the sets 
of n orthogonal vectors in the space having the property that the ith vector of 
each set is a distance +/1/a; from the origin. 





4. Concluding remarks. The importance of the result given in this paper for 
formulas of degree 2 is that it is the first result (other than the trivial one point 
formula, the centroid of R, which integrates any linear function) which holds for 
an arbitrary region in n-dimensional space and which gives all such formulas con- 
taining the minimum number of points. 

A question, which may have some practical importance, which may be asked 
about the above formulas of degree 2 concerns the conditions R must satisfy, say 
for w(x) = 1, in order that such a formula will exist with all of its points interior 
to R. For example, can a formula interior to R be found if R is convex? if R is 
star-like about its centroid? 

The error bound of von Mises [6] for n-dimensional integration formulas is 
very well suited for use with the formulas developed in this paper. In a later paper 
we will give specific values of this error bound for various known formulas. 
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I am especially indebted to Dr. P. C. Hammer for many discussions concerning 
this subject. 
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Calculation of Transient Motion of 
Submerged Cables 


By Thomas S. Walton and Harry Polachek 


Abstract. The system of nonlinear partial differential equations governing the 
transient motion of a cable immersed in a fluid is solved by finite difference 
methods. This problem may be considered a generalization of the classical vibrating 
string problem in the following respects: a) the motion is two dimensional, b) 
large displacements are permitted, c) forces due to the weight of the cable, buoyancy, 
drag and virtual inertia of the medium are included, and d) the properties of the 
cable need not be uniform. The numerical solution of this system of equations 
presents a number of interesting mathematical problems related to: a) the nonlinear 
nature of the equations, b) the determination of a stable numerical procedure, and 
c) the determination of an effective computational method. The solution of this 
problem is of practical significance in the calculation of the transient forces acting 
on mooring and towing lines which are subjected to arbitrarily prescribed motions. 


1. Introduction. This problem arose as a result of an urgent requirement by the 
Navy in connection with a series of nuclear explosion tests which were conducted 
in the Pacific. In preparation for these tests a number of ships were instrumented 
and moored at specified locations from the explosion point. These positions had 
to be maintained intact during the period preceding the explosion. However, the 
bobbing up and down of the ships due to ocean waves could excite transient forces 
in the mooring lines sufficient to break them and thus result in the loss of informa- 
tion from the tests. Several months prior to these tests a request was made to the 
Applied Mathematics Laboratory to calculate the magnitude of the forces acting 
on the mooring lines for waves of varying amplitude and frequency. The two factors 
which made a theoretical solution feasible at this time, whereas it would not have 
been possible several years ago, were: a) the availability of a high-speed computer 
and b) the recent progress made in the understanding and development of nu- 
merical methods for the solution of systems of partial differential equations by 
finite-difference methods. 

Although this problem was solved to satisfy a specific request, it is more useful 
to regard it as the general problem of the two-dimensional motion of a cable or 
rope immersed in a fluid, and it becomes immediately apparent that its solution is 
applicable to a wide class of engineering problems involving the motion of cables, 
such as: a) the laying of submarine telegraph cables, b) the towing of a ship or 
other object in water, or c) the snapping of power lines as a result of transient 
forces caused by storms. The problem may be stated abstractly as follows: Given 
the initial conditions (i.e., position and velocity at any time, f&) and boundary 
conditions (positions of end points at all times) of a cable immersed in a fluid, 
determine its subsequent motions. The motions are assumed to take place in two 
dimensions. 


Received September 24, 1959. This work was done at the Applied Mathematics Labora- 
tory, David Taylor Model Basin, and the paper is a condensation of DTMB Report 1279 [1]. 
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Forces that are assumed acting on the cable are: a) foreed motion of the ex- 
tremities of the cable, b) damping or drag as it moves through the fluid, c) inertial 
reaction of the surrounding fluid, d) weight of the cable, and e) buoyancy. Vari- 
ations in the mass as well as other physical properties of the cable along its length 
are allowed. However, in the present solution it is assumed that the cable is in- 
extensible. The displacements may be large and the motions rapid, provided that 
all significant components of the driving motion lie in a frequency range well 
below the lowest natural frequency of the line for elastic (longitudinal) vibrations. 
In other words, the cable must be sufficiently short (or the velocity of propagation 
of elastic waves sufficiently great) that the line may be considered to be in equilib- 
rium as far as longitudinal waves are concerned. In subsequent work the authors 
have carried out solutions for cables with elastic properties. 


2. Derivation of Equations of Motion. The problem under consideration is a 
generalization of the classical problem of the motion of a vibrating string. We wish 
to deduce the approximate motion of a flexible steel cable without becoming in- 
volved in the explicit computation of the elastic forces which act on the cable. 
However, the formulation of the problem will be more general in several respects, 
namely: 

a) Longitudinal as well as transverse motions of the line must be taken into 
consideration. 

b) The occurrence of large displacements from the equilibrium configuration 
of the line must be permitted. 

c) The extremities of the cable may be at different levels with the cable sagging 
between the positions of support. This requires that the weight of the line be taken 
into account. Thus, even when the line is in static equilibrium, the tension will not 
be uniform nor will the line be straight. ° 

d) Since the cable is submerged, the static forces must include the buoyancy 
of the medium, and the dynamic forces must allow for the virtual inertia of the 
medium. Furthermore, it is necessary to make provision for damping forees due 
to the drag on the line whenever transverse motion is occurring. 

e) Finally, it is desired to suspend concentrated loads at one or more points 
along the line and to change the diameter and linear density of the cable at specified 
points. 

The best approach to the solution of a problem with such general specifications 
appears to be a numerical method based on finite-difference approximations. 
Inasmuch as we are committing ourselves to the eventual! use of differences in both 
the time and space dimensions, it will be simpler to introduce the spacewise dis- 
creteness into the original formulation of the problem. We therefore proceed at 
once to the derivation of the equations of motion of a simplified model in which the 
distributed mass of the cable has been replaced by a series of discrete masses m; 
attached to a weightless, inextensible line. This leads to a system of ordinary differ- 
ential equations. It may be shown that, in the limit, the resulting equations pass 
over into the corresponding partial differential equations for the motion of a sub- 
merged cable. 

Figure 1 shows a typical configuration of the system with the cable attached to 
a float at the surface and anchored to the bottom. Also, a heavy load is suspended 
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Fic. 1. —Mooring line. Fig. 2. —Discrete representation. 


from a point near one end of the cable. Other boundary conditions are possible, but 
the equations of motion will be the same in any case. The horizontal and vertical 
coordinates of a point on the line are called z and y, respectively, and the angle 
between the horizontal and the tangent to the line is designated by @. Figure 2 
illustrates the corresponding discrete model for which the equations will actually 
be derived. The line is divided into segments in such a way that there will always 
be an integral number of them between any points where an abrupt change in some 
parameter occurs. The junctions between the segments are numbered according to 
the subscript index j, which runs from 0 at the anchor to s at the surface. 

Before we can properly invoke Newton’s law of motion, it is necessary to con- 
sider the inertial properties of the fluid in which the cable is immersed. We shall 
assume that the kinetic energy imparted to the surrounding medium is independent 
of the component of velocity parallel to the line, whereas it varies as the square of 
the component of velocity at right angles to the line. Thus, when an element of 
the cable is accelerated longitudinally, no hydrodynamic reaction occurs, but when 
the cable is accelerated transversely, it behaves as though it possessed additional 
inertia. For convenience in formulating the equations of motion in Cartesian 
coordinates, the transverse component of acceleration, namely, —# sin 6 + 9 cos 6, 
can be further resolved into horizontal and vertical components (to which the 


corresponding components of the accompanying inertial reaction will be propor- 
tional). That is, 


Horizontal component = # sin” @ — g sin @ cos 6 
yr . oe 2 a2 
Vertical component = yj cos 6 — £sin 6 cos 6 


Thus, each component of the hydrodynamic reaction depends on both components 
of acceleration. In general, the reaction force is not parallel to the acceleration 
vector (except when the tangential component is zero), so that it is necessary to 
regard the inertial parameters of the system as tensors rather than simple scalars. 

The differential equations governing the motion of the jth station on the line 
(see Fig. 2) can be written in matrix notation as follows: 


I; —K;\ (4; F;* 
(2.1) ' i xg 
—K; J; Yj F; 
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where 
I, = m; + 4(e544 sin’ 054; + €;4 sin’ 0:4) + m;* 
Ji = m; + 3(€544 COS’ 0j44 + €;-4 cos 634) + m;* 
K; = 3(¢;+4 sin 0;44 cos 0;4; + e;4 sin 6;4 cos 6;4) 
andt 
m; = dluisaliog + wjals) 
C544 = pk jbisio iss 
m;* = m;* + pV* 
m;* = m* + evs" 
COS O54; = (Xj41 — 25)/Liay 


sin 0345 = (yin — ys) /liay- 


Each lumped mass, m; , has been expressed as the average mass of the two seg- 
ments of cable which lie on either side of station 7. Also, one-half the equivalent 
transverse mass, e;, of the fluid entrained with each of these segments has been 
included in the inertia tensor. Furthermore, at those stations from which a weight 
is suspended, the effective horizontal and vertical masses, m;* and m;", of the weight 
are to be added. For simplicity in allowing for virtual inertia, we have assumed 
that any such weights possess a certain degree of symmetry and remain upright 
as the line moves about. 

The force vector, F , on the right side of eq. (2.1) can be expressed as the sum 
of internal forces (the tensions acting between adjacent mass elements) and what- 
ever external forces are present. Thus, in expanded form the equations of motion 
can be written 


(2.2) I jé; — Kyjj = T5445 cos 0544, — Tj4 cos 054 + X; 


— Kj; + Jijj = Ti+ sin 034, — Tj4 sin 054 + Y; 


where 7';,; = tension in segment of line between stations 7 and j + 1 
X; = horizontal component of resultant external force at station j 


Y; = vertical component of resultant external force at station 7. 


There are two sources of external force, namely: 1) gravity, which gives rise to the 
weight minus the buoyancy and acts only in the vertical, and 2) fluid resistance, 
which gives rise to the damping forces. Thus, we write 


X; = —3[Dj+ sin 654, + Dj sin 0;4] + X;* 
3[D j44 cos 0;44 + Dj cos 6;4] + Y;* — W; — W;* 


(2.3) 


j 
where 
W; = mig — 3g(ljsyoiag + Laos) 
W ;* = m;*g — pgV;* 


t See Appendix for definitions of the parameters which appear in the formulas. 
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and D;,, = drag on segment of line between stations j and j + 1 
X ;* = horizontal component of damping force on weight at station j 
Y ;* = vertical component of damping force on weight at station j. 


Again, in order to get the best approximation to the continuous case, the net effect 
of the drag at station 7 has been expressed as one-half of each component of the 
drag on the segments which lie on either side of this station. The buoyant force of 
the displaced fluid has been treated likewise. 

We have assumed that the drag, D;,; , on a segment of the line acts in a direction 
at right angles to the line. This is a good approximation whenever the velocity is 
high enough to produce significant forces, since at all but the lowest Reynolds 
numbers the tangential component of the hydrodynamic force is very small com- 
pared to the normal component. Furthermore, we assume that the drag is pro- 
portional to the square of the component of relative velocity normal to the line: 


(2.4) Diss = —Sisa9s44 | a4 | 
where 
Sis = B0CPaliadixs 
Ging = —31 (Zin. — c) + (4; — c)] sin 9543 + 3[¥sar + 9] cos 9544 - 


The positive normal to the line has been arbitrarily taken to be directed upward 
when @ equals zero. The use of the minus sign and the introduction of the absolute 
value of one of the velocity factors ensures that the drag will always be opposed 
to the direction of q;,; and thus act as a dissipative force to remove energy from 
the system. Since the velocities of the two endpoints of each segment will, in general, 
differ slightly, their mean value (which for a straight line segment is exactly equal 
to the velocity of the midpoint) is taken as a representative value in the definition 
of g;+; . In addition, the definition allows for the presence of a uniform horizontal 
current, c, to incorporate the ability to treat towing lines as well as mooring lines 
(or mooring lines subjected to ocean currents). 

In addition to the drag on the line itself, there will also be resistance to the 
motion of any concentrated loads which may be suspended from the line. These 
additional damping forces will vary with the velocity but will not, in general, be 
directed exactly opposite to the motion of each weight. However, on account of the 
assumed orientation and symmetry of any such weights, the resistance force will 
be parallel to the velocity vector whenever the relative motion is either purely 


horizontal or purely vertical. Accordingly, the two components of resistance may 
be written 


. Xj#* = —fj*uj(z; — ¢) 
- Y* = —fj'ug; 
where 

f° = 400;*S;* 
fi" = 40C;"8;" 
uj = [(#; — ¢)’ + 9} 
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Up to this point an explicit formula has been given for the evaluation of every 
term in the equations of motion (2.2) with the exception of the tensions. To de- 
termine these we must invoke the inextensibility condition which was assumed 
at the outset. This takes the form of a constraint on the motion of the line. It re- 
quires that the separation between adjacent stations must not change with time. 
Thus, we write 


(2.6) (x; — a2)” + (y; — ys)? = Gy = const. 


This holds for each segment of the line, and we require that the corresponding 
set of tensions, 7'’;_; , take on values such that the resulting solution of the equations 
of motion will be consistent with eq. (2.6). Because of the implicit nature of this 
condition, we are led to a system of algebraic equations for the determination of 
the proper tensions. At the extremities of the line (j = 0 andj = s) x; and y; must 
be obtained from the boundary conditions, namely: 


mH = Xo(t), Yo = yo(t) 
(2.7) 

x, = 2,(t), Ys = Y(t). 
These are given as functions of time, and permit the introduction of any desired 
types of driving motions. 

Finally, to complete the formulation of the problem a set of initial conditions 
must be given for each station on the line. Since the equations of motion are of 
the second order, it is necessary to specify both the coordinates and the velocities 
at t = 0. That is, 


(0) =2;, y0) = y; (j = 1,2,---s—1) 


(2.8) 
#(0) =a), 9;(0) = 9; (j = 1,2,-+-s—1) 


where the superscript index “0” is used to designate a value at the origin in time. 


3. Solution of Equations by Finite Differences. 


A. GENERAL DESCRIPTION OF COMPUTATIONAL PROCEDURE. The equations 
governing the motion of a cable, as derived in the last section, are summarized 
here. The basic equations of motion, (2.2), are repeated for convenience, 


I #; — Kj; = 1544 cos 054, — Tj cos 0:4 + X; 


(2.2) ‘ 
—K jk; + J ij; > T 544 sin 544 - Td sin 054 + Y; 


2 


where 
a) sis the number of segments into which the cable is divided, 
b) I;, K;, J; are given in eq. (2.1) and are functions of the physical properties 
of the cable and of position only, 
c) X;, Y; are given by eqs. (2.3), (2.4), and (2.5) and are functions of the 
physical properties of the cable and of position and velocity. 
In addition, the motion is governed by the condition of inextensibility of the cable, 
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eq. (2.6), namely, 


(2.6) (x; — aja)? + (yj — yin)? = G4 = const. (j = 1,2, --- 8). 
The differentiated (with respect to time) forms of this relation, 
(3.1) (25 — Xj) (Xj — Zp) + (y5 — ys) (95 — Fin) = 0 

(j = 1,2, 8), 


(3.2) (tj — Xja) (45 — Ea) + (ys — Yi) (G5 — Gia) 


+ (4; — Za)? + (9; -—ga)y=0 (Gj =1,2,---8) 
are also used in the computation. 


For numerical solution by finite-difference methods the following finite-difference 
equivalents are used, 





n+l n n+l n 
= n+} x — -nt+} y a : ‘ 
(3.3) so = At ms lees coms” (j = 1,2,---s—1), 
etl — +. oi n+l 9 n n—l 
a” os vj 60; + aj ij;" = Yi 245 TF Yi 
(34) ° (At)? om (at)? 


(j = 1,2,---8s—1). 


It is assumed that the boundary and initial conditions are known. These are given 
in eqs. (2.7) and (2.8), respectively. The system of equations summarized above, 
consisting of eqs. (2.2), (2.6), (2.7), (2.8), (3.1), (3.2), (3.3), (3.4) with the 
auxiliary equations (2.1), (2.3), (2.4), (2.5), completely describe the motion of 
the cable. 

The computational procedure, as developed in detail in the remainder of this 
section, consists of an algorithm to determine the values 2}*’, y?*', #7, g?™ 
(at time t = t"*’ = t" + At and t"* = ¢” + 44t) from known values z;”, Yi", 
é? +, yf? (at time t = t" and t””). It is convenient to divide this algorithm in 
two phases. The first phase involves the determination of a tentative (but con- 
sistent) set of tensions for all segments, and numerical integration of the equations 
of motion to predict the coordinates one step ahead; the second phase involves the 
evaluation of small discrepancies in the constraint equations from which a set of 
first-order corrections to the tensions can be obtained, and integration of the equa- 
tions a second time to obtain accurate values of the coordinates. 

Phase 1. Using eqs. (2.2) and (3.2), (3s — 2 equations), we compute the (3s — 


2) 
unknown variables #;", j;" (j = 1, 2,---s — 1) and T}-y (j = 1, 2,--- 8). 
We now use eqs. (3.3) and (3.4), (4s — 4 equations), to compute the (4s — 4) 
variables at the next time step, 27’, y?', 27", gf? (j = 1,2 --- s — 1). These 


are considered to be only tentative values (denoted in subsequent text by use of 
the tilde). 

Phase 2. To obtain improved values of the tensions Tj, (j = 1, 2, --- 8) 
and the quantities #;", g;", 27", y7“, 274, g7* (a total of (7s — 6) quantities) 
we use the system of equations (2.2), (2.6) and (3.3), (3.4), consisting of (7s — 6) 
equations. However, since eqs. (2.6) are not linear but quadratic in the unknowns 
x; and y;, an explicit solution is impractical to obtain. For t»is reason a computa- 
tion algorithm based on the Newton-Raphson method of successive approximations 








34 THOMAS S. WALTON AND HARRY POLACHEK 
is developed. A detailed discussion of the computational procedure used in this 
problem is given in the sections which follow. 


B. DETERMINATION OF TENTATIVE VALUES OF TENSIONS. The system of equa- 
tions (2.2) may be regarded as a set of (2s — 2) linear equations in the variables 
%; and ¥; (accelerations) and may be solved directly for these variables. If we 
designate 


L; = (At) 1;/(1jJ; — Kj) 
M; = (At)’ J;/(IjJ; — K/) 
N; = (At)’ K;/(1jJ; — K?), 


then the equations of motion (2.2) can be reduced to: 
#; = [RjT jx, — PjT 34 + Uj)/(At)’ 
(SjT j43 — QiT 34 + Vi/(At)° 


Yi 


P; = M; cos 0;4 + Nj; sin 0; 
Q; = N; cos 6,4 + LZ; sin 0;4 
R; = M; cos 0j4, + N; sin 04, 
S; = N; cos 0;4; + LZ; sin 054; 
U; = M; X; + N; Y; 
V; = N; X; + L; rs. 

We observe that eq. (3.2) involves positions, velocities’, and accelerations. As is 
often the case with finite-difference procedures, it proves to be convenient to com- 
pute positions and accelerations at the mesh points while velocities are found at 
the midpoints in time. For this reason we shall use a modified form obtained by 


evaluating eq. (3.2) at ¢ = ¢” and at t = t””’, and then adding the two results 
together, namely, 


(x5” — 25-1) (2;" — 25) + Cys” — ys) (99° — F5-) 
(36) 9 + Ft = FS) a - 8) + WI a) rt - 
+ 2(At)™ [(2;* — 27") — (aja — 2f5)P 
+ 2(At)™ [(y;" — y7*) — (yf — yFR)P = 0, 
in which we have used the approximations, 
aj" — ofa)? + (a7 — AFT)? we 2(a77 — aR)? 
FY A(x" — a7*)/At — (aja — aja) /Aty, 
(g3" — Ga)? + (OF — GRY? & 20g77 - BAY’ 
 2(y;" — yj ')/At — (yf — yfa)/Aty’. 


Note that eqs. (3.6) are linear in the accelerations. Likewise, eqs. (3.5) are linear 
in the tensions. Consequently, when these latter expressions are substituted for the 
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acceleration components in the constraint equations (3.6), we obtain a set of condi- 
tions which are linear in the tensions, namely, 


E34 T3-+4 — Fi Tha + Gia Thy 
+ ES TH — FSS TH + GS TH + Hi + HS 


(3.7) ” 
+ 2(zj* — 23“) — (af — zfa)P 
+ 2(yi" — x3") — (ya — ya )P = 0 
where 
Ej-4 = (x;" = 2j-1)Ph- + (y;" - yj3—-1) Q5-1 


F5-y = (a;" — ap) (P;* + Rha) + (y7* — yj) (Q;? + Sj) 
G54 = (25" — 25) R;” + (y;” — yj) S;" 
H3-y = (2j° — 25-1) (Ui? — Usa) + (ys? — ypu) (VI? — Vi). 


Now assume that the solution is correct up tot = t”. Then all quantities in (3.7) 
can be evaluated at once except for T3_,, Tj, and T7,, . The tentative values of 
the tensions—signified by the tildes—are determined by the following system of 
equations: 





i —Fos Gos 7 I 0.5 7] T—Ws 
Ets —Fis Gis Tis —Vi's 
(3.8) N25 ; —F?2s5 : Gos Ts i —¥is 
Ets Ft; Gras gar —Viis 
q Et_os —Ftos | oan! q —Vi05 a 

















where 
n a—l n—1 n—l n—1 yn—l yn—l1 n—l n 
Wj = Ej Ti — Pig Tig + Gi Ting + Hi + Aj 


+ 2[(a;" — 237) — (apa — aT )P + 2A(y* — yf?) — (ypu — yf)P. 


In general, we can write (for 7 = 1, 2, --- s) 
(3.9) Ej T3-+4 — F3-4 T34 + G4 Ti + ¥ G4 = 0 


with the conditions: Eo; = Gi, = 0 forall n. Also, Po” , Qo” , Ro” , So” , and P," , 
Q.”, R.” , 8S." = 0 for all n; and 


U," = (At)*é" and U,” = (At)*#," for all n, 
Vo" = (At)*go" and V," = (At)*g," for all n. 


The matrix of coefficients of the system of equations is a triple diagonal one, and 
it can be easily reduced to a single linear equation by elimination. Thus, we solve 
eq. (3.9) for T},, . 








no Fit gm, _ Ein qm _ Vins 
(3.10) — Gry Th Gry j-i Gry 
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Now we express each tension as a linear function of 77.5 (the tension in the first 
segment) as follows: 


(3.11) Ti = O44 Tos + BF +4 

and we arrive at the following recursion formulas for a}4; and 6},,, namely, 
aj4y = (F354 aj — Ej aj4)/Gj 

(3.12) 


Bisy = (F5-4 Bi-y — Ej-4 Bj-y — Y¥i-+)/Gi-+ 

with the conditions 

aos = 1, aos = 0 for all n, 

Bos = 0, Bos = 0 for all n. 
Starting with 7 = 1, we evaluate aj,; and 8,; recursively up to7 = s — 1. We 
then find 7¢.5 from the last equation of the system. 

zn F-4Bs-4 — Et+4 Bry — Vi 

ae) _ desc. a nz F 





C. Computation oF TENTATIVE CoorpiNaTEs. In order to solve eqs. (3.5) 
° ~ 1 ~ 1 . . . . 
numerically, we replace #;' and #; by their simplest central-difference approxima- 
tions, eq. (3.4), namely, 
? #j" = (aj — 2aj" + 27") /(At)" 
(3.14) . _— 


gi" = (yp — Qy;” + yF*)/ (at). 


n 


7 n +1 +1 . ° ° ° 
Now we solve for x; ~ and y;* , considering these as tentative values subject to a 
slight modification in order to satisfy a system of constraints. Thus we write 


n+ ‘ —l Y 
Ep = Qe — af — PP Th + RP Thy + *UF 
an+l 


Gj) = 2y — yf — OF TH + SPT + Vi ~ 


n 


The quantities P;", Q;", R;", S;", U;" and V,;" are the same as were used to set 
up the coefficient matrix for the tensions, and the values for T7_; and T7,; are 
obtained from (3.11). 


(3.15) 


D. DETERMINATION OF IMPROVED VALUES OF TENSIONS. Next, we determine 
the set of corrections 57’; to be applied to the tensions T;, in order that the 
values of z}** and y?** should also satisfy the inextensibility condition (2.6). For 


this purpose we define the function 
(3.16) ort = (ap — cht)? + (yf? — yy -— Gal 


which measures the discrepancy in the distance between the extrapolated positions 
of pairs of adjacent stations. We observe from eqs. (3.15)—with the tildes sup- 
pressed—that z?*' and y?*’ are functions of the tensions. Consequently, 2"“} may 
also be expressed as a function of the tensions. This enables us to write the 
system of constraints to which the tensions are subject as follows: 


(3.17) OFF = OFF (T34, T4,Tit =0 8 (f = 1,2,--- 8) 


since 277} vanishes when the inextensibility condition is obeyed. 











TRANSIENT MOTION OF SUBMERGED CABLES 37 








Now let 
(3.18) Tia = Thy + 8734, 
and expand 27+} in a Taylor series about the point {T7-4, T7-+4, Tru}. Thus, we 
obtain 
" - an7} . ann} 
ort = OF + aT sT744 + ate, éT 74 
(3.19) nt! 
+ a 5T 544 + higher order terms, 
where 


OFF = OFT, Th, Thal 
= HCE" — EPL)’ + GF" — GAY — G4). 
Provided that the tentative values 7}, are sufficiently close to the correct values 
Tj; , we may neglect the higher order terms in the expansion (3.19), and thereby 
obtain a system of s linear equations for the differential corrections 5T}_, . These 


equations have the same form as the previous system (3.8) for determining 77_, , 
namely, 


























[—Pet! Osi" 1 rere. 7 [agit 7 
n+l a” n+l n+l ™ nnad 
1.5 1.5 1.5 6T is —Ors 
n+l n+l n+l ». 
(3.20) 2.5 —— 2.5 ; G23 &T es oat — ie?" 
nay —Frt. otis éTe-15 —OP*!, 
a B82, —Frts | Tos | a —-. J 
the general expression being (for 7 = 1, 2, --- s) 
(3.21) py OT — PRAY OTT, + GU} oT74 + GE = 0 
where 
: an; - ~ - ieee 
Ey = aT, ; = (877 — SPL + (GF -— FADO 
a 
pet e077} antl ze ~n+1 ~n+1 n n 
i = oT : = (Z; —_ 1)(P;* + R}-1) + (9; — Yj-1 )(Q; + Sj-1) 
i- 
n ann} ~n ~n ~n ~n yn 
G7? = aT" _ (z} ee ay) R? + (7; +1 — 9777)S; 
i+ 
ay 


with the conditions: E¢{* = 
S;” being the same as in (3.7). 

The system (3.20) can be solved in a manner completely analogous to the solu- 
tion of the system (3.8). Thus, we write 


= 0 for all n, and the quantities P;*, Q;", R;* and 


(3.22) bT 544 = «3448705 + Ajay 
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and obtain the following recursion formulas: 


4 n+1 pint. ra +1 
ai Kian = (Pap — By h4)/GP 
(3.23) n n+1 rn+1 ptt 1 
n " Q” Ayn + 
Aj4 = (Pry nj-y ss Ej )j , )/G; 
with the conditions 
kos = 1, kos = 0 for all n, 


dos = 0, dos = 0 for all n. 
Finally, the last equation of the system enables us to solve for 675.5. The result is 
(Prt ay — Bray — ont) 
( —. oa Br n 3) 
We can now obtain the corrected values of the tension in each segment. Thus, 
= 
Tj = Thy + 6Tj 
= Thy + «i6T0s + APY 





(3.24) 6To5 = — 


(3.25) 


E. Computation OF ImMpROvED CoorpinaTeEs. The corrected values of the 
coordinates are found using eqs. (3.15)—but this time with the tildes suppressed— 
namely : 


+1 “i T rm 

x ; = 22; — % C; — P; ‘ ty — R;" Ti 44 oe U; " 
n+l ‘ n n—l Yn n rn 
yi = 2y; — fF — QF Ti + S/T iy + Vi 


For solution on an automatic computer, it is more convenient to express eqs. (3.26) 
in terms of corrections to be added to the tentative values of the coordinates. That is, 


(3.26) 


6x} +1 “a me /sTF- 4 R ET 
(3.27) n+l n n n n 
by; = —Q,ST4 + SST j4. 
Then the corrected coordinates are given by: 
ae -_ ast + ba?" 
(3.28) 
ti ‘ais gt + by? 


F. SpecraL Form or Equations For ComputinG First Time Step. We assume 
that the initial velocity components are zero at each station, and we obtain the 
initial coordinates from the equations for static equilibrium of the line. Since z, 
and 4; = 0, eq. (3.2) reduces to 


(3.29) (xp — x$4)(4? — a) + (vy? — yada? — Hf) 


and, on substituting the expressions (3.5), we find that the tensions are subject 
to the constraint 


(3.30) Ej a — Fi4 Ha + G4 Ma + Hi = 0. 


Comparing this with eq. (3.9), we see that 


(3.31) Sa = H3_, Gj -_ ¥ A eee 8)- 
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The system of equations (3.8) is then solved in the usual way to get the proper 
initial tensions T%4_, . 
To obtain tentative values for the coordinates at t = ¢', we make use of their 
Taylor series expansions about the point ¢ = ¢’, namely: 
vj; = xj + (At)e? + $(At)*e? + --- 
yi = yi + (At)y? + 3(At)*G? + --- 
Taking z;’ and y;’ = 0, and substituting eqs. (3.5) for #; and g,’, we find 
Ei = 2; +3 [— PPT 4 + RPT 4, + U7 
~ 1 0 0 1, O70 , 0 
9; =y¥i + 3 (- Q; T4 + 8; T54 +7 j |. 
The corrections to the tensions are then determined by the system of equations 
(3.20) in the usual manner. Finally, the corrections to the coordinates are com- 
puted as follows: 


(3.32) 


(3.33) 


bri = 4([— P} 8T§4 + R? 6T} 44) 
by; -_ 3 [- Q; 6T)4 + S; 6T% 43] 


and the corrected coordinates are given by: 


(3.34) 


Xj = #; + ba; 
(3.35) ends Se ; 
Yi = 9; + dy;. 

4. Analysis of Numerical Stability. In order to obtain a valid solution of the 
system of partial differential equations governing the generalized motion of a 
cable, it is necessary to insure the stability (in the sense discussed in [2], [3], [4]) 
of the equivalent finite-difference system (2.2), (2.6), (3.3), (3.4). In this section 
we will derive the criteria for stability of this system of equations. We will also 
show that whereas the system of finite-difference equations (2.2), (2.6), (3.3), (3.4) 
is stable for sufficiently small time intervals At, the system (2.2), (3.2), (3.3), (3.4) 
is always unstable. This characteristic of the latter system has led to the abandon- 
ment of this simpler set of equations in favor of the more difficult nonlinear system 
(2.2), (2.6), (3.3), (3.4). 

In order to determine the stability of a system of finite-difference equations, we 
study the growth of a small disturbance or perturbation. The conditions for sta- 
bility are said to be satisfied if the amplitude of a small disturbance, introduced 
at any time, ¢, in any of the dependent variables, does not increase exponentially 
with successive time steps. This condition may be stated as follows: 

If 5F(s, ¢) and 6F(s, t + At) are values of a variation (or perturbation) in any 
of the dependent variables x, y, T in the system, then it is said to be stable provided 
| 6F(s, ¢ + At)/éF(s, t) | < 1. We introduce perturbations 62x, dy, 5T in the 
dependent variables x, y, T, respectively. To simplify the stability investigation, 
we shall omit the terms pertaining to the suspended weights as well as all terms 
involving virtual inertia. From eqs. (2.2), (2.6), (3.3), and (3.4) we obtain the 
following variational system of equations: 


mz; = T 544 6 cos 0544 _ T 54 6 cos 654 + cos 0544 5T 544 
— cos 6;4 6T 4 — 3[Dj4, 6 sin 6;4, + Dj, 6 sin 6;_, 
+ sin 654; 6Dj4, + sin 0;4 6D;-4), 
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(4.1) my; = T 54,6 sin 654; — Tj4 6 sin 054 + sin 0j4; 67 544 

— sin 0;~; 67 ;_; + 3[Dj4, 6 cos 6,4; + Dj+6 cos 6;_, 

+ cos 054; 6Dj4, + cos 6;_; 6D;}, 

cos 654; 6 cos 054; + sin 054; 6 sin 0j4; = 0, 
where 
SDj44 = —2firs | ai44 | Oa» 

qi43 = —31(Zinn — €) + (4; — €)] 8 Sin 0545 + 2°94. + G5) 6 COS 8545 

— ¥ Sin 9543(6%j41 + 64;) + 3 COS O545 (Hj. + 595); 
and where 

§ cos O54, = (8x41 — 6x;)/bi4g, 8 Sim Oj4y = (Syjar — bys) /Ls44 5 
oa7* = (d;" — bef *)/At, dy * = (dy;" — dy7*)/At; 
ba," = (dxf 7" — 26x;" + dx7*)/(At)*, yj” = (SyF** — dy,” + dy") /(At)*. 
We will assume in this analysis that within a small region in the (s, ¢) plane the 

coefficients (7';", cos 6;", D;”, etc.) of the variational functions vary only slightly 
and hence may be treated as constants. We will denote these simply by 7’, cos 86, 


D, ete., omitting the indices. A solution of the system of equations (4.1) can then 


be obtained in the form oe 
éx;” ~~ etitanas 


by;" - b ePitandt 


j At 
8T;" meg eibitan 


where a, b, c are real constants and a complex. Substituting in eq. (4.1), we obtain 
a system of linear homogeneous equations for the quantities a, b and c which has a 
non-trival solution, provided the determinant of the coefficients is identically zero. 
After some algebraic simplifications, the determinant of the coefficients may be 
written in the form 


| F—Asin@ D’+Bsin@ cos@| 
(4.2) —D’+Acos@ F — Beosé sin 0) = 0 


cos 6 sin @ 0 
where 
A = f\|q| (2iy sin 8 — (I sin 6/At)(1 + cos B)(1 — X")] 
B = f\q| (2i(z — c) sin 8 — (1 cos 6/At)(1 + cos B)(1 — X”*)] 
D’ = iD sin 8 
F = mlt/(At)’ + 4T sin’(8/2) 


and where 


A=e™' gE =(A— 24d"). 
Multiplying the elements of the determinant and simplifying, we obtain 
(4.3) A sin@+ Beos#@-— F=0 
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But, 


A sin 6+ Bcos 6 = f|q|[(2isin 8B) p — (l/At)(1 + cos 8)(1 — X")] 
where 
p = ( — c) cos 0+ ysin 0 
i.e., the tangential component of the velocity of the cable (relative to the medium). 


Substituting in equation (4.3), we finally obtain the characteristic equation of the 
variational system (4.1), namely, 


mir’ + {f | q| At (1 + cos 8) — pAt(2i sin 8)] 
+ 4T(At)* sin*(8/2) — 2ml}x + [ml — f | q| lAt(1 + cos B)] = 0. 
Now, comparing the first and second terms of the coefficient of \, we find that the 
second term is negligibly small provided 2pAt < l, i.e., the tangential distance 
traversed by the cable in one time step is very small compared with the length of 
the cable segment. Since this is usually the case and, at any rate, can always be 


satisfied by taking the time step sufficiently small, we will omit this term from our 
subsequent analysis. 


For the case of negligible drag, i.e., f = 0, approximately, we obtain from eq. (4.4) 
(4.5) »? + [47 sin® (8/2)(At)’/ml — 2] + 1 = 0. 


In order for the solution to be stable, the conditions | 4, | < 1 and |A2| < 1 must 
both be satisfied. But if \; is a solution of (4.5), then A» = 1/A, is also a solution. 


It follows that the conditions for stability can be satisfied only if | Ax | = | 1/A, | = 
|r: | = 1. Now, let A; = cosy +isiny, Ax = cos y — i sin y = 1/A,; then, 
| Ar + A2| = | 2 cosy| S 2. Furthermore, from eq. (4.5) we have 


_ 4T sin*(8/2)(At)* 


Ai + he = 2 
ml 





We thus obtain the inequality 


|2 — 4T sin’(B/2)(At)’/ml | < 2 


’ 


or 
(4.6) T sin’(8/2)(At)*/ml < 1. 


This requirement is tantamount to the condition 





ats 4/m@ = 
7 T velocity of transverse wave 


In the more general case, allowing for finite drag, eq. (4.4) may be reduced to 
(after neglecting the second term of the coefficient of \) 


mln’ + [2f | q | Lat cos” (8/2) + 4T (At)? sin’(8/2) — 2mljr 
+ {ml — 2f | q| lAtcos® (8/2)] = 0. 


This equation is more difficult to analyze. However, it is possible to show that both 
|¥:| S 1 and] A:| S 1, provided that we have chosen 
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ml 


Ats a 
= T? 


and 
(4.8) 
m 

f\a\- 
These conditions (4.8) are both necessary and sufficient to insure stability.t 

We will now show that the replacement of eq. (2.6) by its differentiated form 
(3.2) results in an unstable system; and that, furthermore, the use of any time 
interval At, no matter how small, does not change the unstable character of the 
equations. It will suffice to show that this condition exists in the case when the 
drag is negligible, i.e., f = 0. The variational equation corresponding to eq. (3.2) is, 


At < 


(x, — 2j1)(8%; — 6454) + (4; — £3) (6x; — S254) 
+ 2(4; — 25-1) (84; — 64541) + (ys — Yin) (69; — Hj) 
+ (95 — Gi-r) (Oy; — Syjn) + 2095 — Yin) (89; — SYju) = O. 


Substituting appropriate values for 6x and éy and neglecting terms containing f, 
we find that the determinant equation (4.2) is replaced by 





F 0 cos 6 
(49) 0 | F sin 6 0. 
(a; — ajalé + (&; — #4)(At)? (yj — yjpwdé + Gs — Gia) (AD 0 
+ 2(¢; — éja)(1 — d’)At + 2(y; — yia)(1 — “At 


Multiplying the elements of the determinant, we obtain 

F cos O[(a; — aj1)E + (¥#; — €)1)( At)? + 2(4;-— 2;4)(1 — XN") AL 

+ F sin (yj — yia)& + (Gi ~ Gia) ( At)’ + 2(G; — Gia) (1 — N*) AL] = 0. 
Equating 


(4.10 


cos 0 = (2; — 2;-1)/l, sin 6 = (y; — yj-)/1; 
and using the relations (3.1) and (3.2), we obtain in place of eq. (4.10) 
(4.11) FiCE — [(; — #54)’ + (9s — Gea)*"M(At)4 = 0. 


In this case the characteristic equation is a quartic, and it is necessary that the 
absolute values of all four roots be <1 to insure the numerical stability of the 
system. Thus, we must examine the roots of both factors of eq. (4.11), namely, 


(4.12) F=0 
and 
(4.13) Pe — [(¢) — 4)a)’ + (95 — ¥ia)*(At)* = 0. 


It can be shown that eq. (4.12) is equivalent to the criterion (4.6) which can 





+ When the terms for the virtual inertia and suspended weights are included in the stabil- 
ity analysis, it is found that the quantity m in eq. (4.8) should be replaced by the expression 
m + m* + e+ p(V*sin’é + VY cos’). 
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be satisfied provided 


ats ,/m 
a vi 


However, the pair of roots associated with the second factor of (4.11) cannot satisfy 
the stability condition for any finite At because eq. (4.13) requires that 


E = [(@; — j4)° + (G5 — Gea)" at)’/P, 
an inherently positive quantity. This conclusion follows as a result of the definition 
&= ) — 2+". If); is a root of equation (4.13), then 4 = 1/As is also a root of 
this equation. As before it follows that for stability | As | < l and | Ay | = | 1/As| S 1. 
Hence, | As | = | Ay | = 1. Let \3 = cosy + isin y, 4s = cosy — isiny = 1/);3; 
then £ = 2(cosy — 1), or —4 S & <= O. Thus, to satisfy the stability requirement 
— must lie between 0 and —4, and consequently can never be positive. 


5. Numerical Results. A number of solutions were carried out with the anchor 
end of the line fixed and with the surface end forced to follow the motion of trochoidal 
waves of varying amplitudes and periods. Several typical solutions are reproduced 
here for the information of the reader. In Figures 3 and 4, plots are given of the 
maximum tension attained along the cable as a function of time for wave heights 
of 6 feet and periods of 12.5 seconds and 5 seconds, respectively. The periods of the 
variation in maximum tension correspond to the periods of the forced vibration, 
as expected. The maximum tension, however, increases in magnitude from 33,250 
Ibs in the case of waves of 12.5 sec period to 49,500 lbs when the period is 5 seconds. 
In Figure 5, the maximum tension attained for wave heights of 9 ft and a period of 
7.5 seconds is plotted. The maximum tension is approximately 60,000 lbs as com- 
pared with 38,500 lbs for the case of 6-ft waves with the same period. 
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Fig. 3.—Mooring line oscillations: wave height, 6 feet; period, 12.5 seconds. 
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Fig. 4.—Mooring line oscillations: wave height, 6 feet; period, 5 seconds. 
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Fic. 5.—Mooring line oscillations: wave height, 9 feet; period, 7.5 seconds. 


As an experiment to aid in understanding the effect of the drag on the motion 
of the cable, one case was carried out with zero drag (i.e., motion in vacuum). 
A very interesting motion pattern was obtained which does not possess a periodic 
character. The resulting maximum tension is plotted in Figure 6. 

The programming of the various phases of this problem was carried out by Mr. 
Thomas McFee, of the Applied Mathematics Laboratory, in a most effective 
manner. The speed and accuracy with which he accomplished this phase of the 
solution were largely responsible for the success in meeting the required time 
schedules. The authors would also like to express their gratitude to Mr. R. T. 
McGoldrick, of the Structural Mechanics Laboratory, for proposing this problem 
and for a number of helpful discussions; to Dr. R. Bart, Structural Mechanics 
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| Fic. 6.—Mooring line oscillations without drag: wave height, 6 feet; period, 7.5 seconds. 
| 
Laboratory, for a number of ideas used in setting up the numerical procedure; to 
Dr. E. H. Kennard, David Taylor Model Basin, and Dr. R. M. Langer, Bureau of 
Ships, for helpful discussions in connection with the definition of the problem; to 
Dr. Daniel Shanks, Applied Mathematics Laboratory, for valuable suggestions; 
and to Miss Corinne Lundgren, Applied Mathematics Laboratory, for assistance 
in the preparation of the figures. 
Appendix—Notation. 
C?., Drag coefficient for segment of cable between stations j and j + 1 
C;* Resistance coefficient for horizontal motion of suspended prism 
C;" Resistance coefficient for vertical motion of suspended prism 
c Velocity of uniform horizontal current 
50 D Drag 
d;,; Diameter of segment of cable between stations j and j + 1 
€;4; Virtual mass of entrained fluid between stations 7 and 7 + 1 
hie F Resultant force : 
m). Fis Drag factor for cable = (0/2) Cissbiridins By de 
dic fi Horizontal drag factor for suspended prism = (p/2)C; S; 
f;" Vertical drag factor for suspended prism = (p/2)C;"S;" 
Mr. g Acceleration due to gravity 
tive I; Component of inertia tensor 
the i Imaginary unit 
isan J; Component of inertia tensor 
Tv. i Subscript denoting station number along line 
lem K; Component of inertia tensor : 
nies kj, Virtual inertia coefficient for segment of cable between stations j and j + | 
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lj4; Length of line between stations j and j + 1 

m; Mean mass of segments of cable adjoining station j 

m;* Mass of prism suspended from station 7 

m;* __ Effective horizontal mass of suspended prism 

m;” _ Effective vertical mass of suspended prism 

n Superscript denoting time-step number 

0 Superscript denoting initial state (origin in time), or subscript denoting 
anchor end of line 

p Tangential component of velocity of cable (relative to medium) 

q Normal component of velocity of cable (relative to medium) 

S** Projected area of suspended prism along z-axis 

S 


;* Projected area of suspended prism along y-axis 
8 Subscript denoting surface end of line 
T Tension 
t Time 
At Time-step interval 
u Magnitude of velocity of cable (relative to medium) 


V;* Volume of prism suspended from station j 

V;* Equivalent volume of horizontal virtual mass of suspended prism 
V;” Equivalent volume of vertical virtual mass of suspended prism 
W; Mean net weight of segments of cable adjoining station j 


W,* Net weight of prism suspended from station 7 

X Horizontal component of resultant external force 

X,;* Horizontal component of damping force on suspended prism 
2 Horizontal coordinate of cable 

Y Vertical component of resultant external force 

Y;* Vertical component of damping force on suspended prism 
y Vertical coordinate of cable 

a Damping coefficient of the perturbation functions 

B Angular wave number of the perturbation functions 

6 The variation of 

6 Angle between horizontal and tangent to cable 

r Eigenvalue (root of characteristic equation ) 


uj+; Linear density of segment of cable between stations j and j + 1 

p Density of fluid medium 

a;+; Cross-section area of segment of cable between stations j and j + 1 
Dot signifies differentiation with respect to time 

~ Tilde signifies tentative value of a variable 


Applied Mathematics Laboratory 
David Taylor Model Basin 
Washington 7, District of Columbia 
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Abscissas and Weights for Lobatto Quadrature 
of High Order 


By Philip Rabinowitz 


In recent years, Gaussian quadrature has become the standard method for 
numerical integration in many computer installations [1, 2]. In general, Gaussian 
rules are most economical since an n-point rule is exact for polynomials up to 
degree 2n — 1 and no rule can do better. However, for particular classes of func- 
tions and for particular applications, other rules may be more efficient. Thus there 
are occasions when we prefer a closed rule, i.e. one which includes among its ab- 
scissas the two end points of the integration interval. This is the case when the 
integrand vanishes at the two end points as it does in Longman’s method for 
evaluating integrals of oscillating functions [3]. If we want to check a quadrature 
over a given interval by doing two additional quadratures, each over half the 
interval, then a closed rule will save at least two evaluations of the integrand. 

Let us normalize our integration interval to (—1, 1) and consider the closed 
symmetric n-point integration rule with n odd, n = 2m + 1, 


1 m 
(1) [f2) de = Daun fain) + En 
1 —=—m 
with %imn = +1, tin = —2Xen, and Gin = A. It was first shown by Lobatto* 


[4] (cf [5]) that the 2n — 2 quantities, n — 2 abscissas z;,, and n weights a,,, , can 
be so chosen that the integration is exact for polynomials up to degree 2n — 3. 
The abscissas z;, are the zeros of P,_, (x), the first derivative of the Legendre 
polynomial of degree n — 1, P,_:(x), and the corresponding weights a;, are given 
by 
2 

Aken 2 Ss . ar ee 

. n(n — 1)[Pa-i(2en)}? 
while dan = @emn = «/n(n — 1). In this case 


n(n — 1)°2"""*[(n — 2) I° con» 
~@n—Dian—ajp % 


(2) 


Tables of abscissas and weights for Lobatto quadrature (where only the values 
of k = 0,---, m are listed) were calculated on WEIZAC to 19D for n = 5(4) 


Received May 20, 1959; in revised form August 7, 1959. 

* According to Scarborough [6], Lobatto modified Gauss’s integration rule so as to include 
the end values and also the value of the function at the midpoint of the interval. Radau [7] 
extended Lobatto’s method to the general case of m fixed abscissas in an n-point rule (m < n). 
Kopal [8] designates by Radau quadrature all rules with m preassigned abscissas. However, 
Hildebrand [5] and others use Radau quadrature to denote the case where only the left-hand 
endpoint of the integration interval is preassigned and Lobatto quadrature, the case where 
both endpoints are preassigned. Radau [7] tabulated Lobatto abscissas and weights for n = 
2(1)8 to eight decimal places and for n = 9(1)11 to ten decimal places and Kopal [8] has re- 
produced them. 
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25(8) 49(16) 97. Table 1 lists these values, with the exception of n = 81, 97. 
Values for n = 81, 97 have been deposited in the Unpublished Mathematical 
Tables File. The zeros of P’._,(x) were computed with double precision routines 
using Newton’s method, 


TaBLe 1 





Abscissas and weights for Lobatto quadrature 








akn 





n=5 
-0000000000000000000 -TUIMILILIL 111111 
-6546536707079771438 .5444444444444444444 
1.0000000000000000000 -. 1000000000000000000 
n=9 
-0000000000000000000 .3715192743764172336 
.3631174638261781587 . 3464285 109730463451 
.6771862795107377534 . 2745387 125001617353 
. 89975799541 14601573 . 16549536 15608055250 
1.0000000000000000000 -.0277777777777777778 
n= 13 
-0000000000000000000 . 2519308493334467360 
. 249286930 1062399926 .2440157903066763565 
-4829098210913362017 . 2207677935661 100861 
.6861884690817574261 . 1836468652035500920 
-8463475646518723169 . 1349819266896083491 
-9533098466421639119 .0778016867468189278 
1.0000000000000000000 .0128205128205128205 
n=17 
-0000000000000000000 . 1906618747534694333 
.1895119735183173883 . 1872163396776192359 
.3721744335654770419 . 1770042535156578704 
.5413853993301015391 . 1603946619976215395 
.6910289806276847054 . 1379877462019265591 
.8156962512217703071 . 1105929090070281614 
- 9108799959 155735956 .079498270503687 1192 
-9731321766314183142 .0449219405432542096 
1.0000000000000000000 .007352941 1764705882 
n = 21 
.0000000000000000000 . 1533851903321749485 
. 1527855158021854660 . 1515875751116813845 
.3019898565087648873 . 1462368624479774593 
.4441157832790021012 . 1374584628600413436 
. 57583196026 18306869 . 1254581211908689480 
-6940510260622232326 . 1105170832191233353 
. 7960019260777 124047 .0929854679578860653 
.8792947553235904644 .0732739181850741442 
, .9419762969597455343 .0518431690008496251 
.9825722966045480282 .0291848400985054586 
1.0000000000000000000 .0047619047619047619 
n = 25 
-0000000000000000000 . 1283083892986619283 
. 127957059483 1069727 . 1272549775383314470 
. 2538130641688765802 . 1241120389379502907 
.3755014578592272332 .1189311794068118254 
. 4910241 148188783826 .1117974662683208882 
.5984841472799932681 . 1028280303479578308 
-6961170488151343668 .0921701399106204219 
. 7823196592407 167804 .0799987748362929818 
.8556764658353165775 .0665137286753127847 
.9149827707346225783 .05193622836849 14746 
.9592641382525344789 .03650473879427 13720 
. 98778994493 14937093 .0204651689329743853 
1.0000000000000000000 .0033333333333333333 
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-0000000000000000000 
.0965481881761070063 
-1921949314674772257 
- 2860472014876740416 
-3772287242533936350 
-4648881616321067560 
-5482070599191116231 
-6264074912812682573 
-6987593166181625957 
- 7645870017935286280 
-8232759230040674696 
-8742781007505622206 
-9171173034509412408 
-9513934513969957434 
-9767861633169063015 
.9930563584336583437 

0000000000000000000 


.0000000000000000000 
-0775101379362542756 
. 1545541206838400005 
. 2306685965595302588 
.3053958040294457421 
- 3782863247310832047 
-4489017863148284990 
. 516817498847 1293170 
- 5816250089 152126880 
. 6429345560645854229 
. 7003774167806829578 
. 7536081218869513649 
. 80330653395809 17814 
.8461797721039918422 
-8849639721696151228 
-9184258706979038132 
-9463641996474921212 
-9686108695676818301 
-9850318591247247763 
-995527 1274452843615 

0000000000000000000 


-0000000000000000000 
.064740137990899 1291 
.1292087333163799330 
. 1931353822540207 100 
- 2562519541901461885 
. 3182937 162511015596 
. 3790004436827688739 
-4381175113207418406 
-49539696 15727524004 
-5505985444332280122 
603490725 1667665006 
-6538516554322124686 
.7014701037710529541 
-7461463415517796826 
. 787692980643654009 1 
.8259357592850707376 
. 8607 142728047894851 


891882646026442427 1 


-9193101442725844907 
-9428817196631997580 
-9624984879932527299 
-9780781245013720159 
-9895551246192307357 
- 9968804559 100229277 
1 .0000000000000000000 


n 


33 


49 





49 





akn 


-0966987 101027 165590 
-0962472849729854620 
-0948972243945918158 
-0926611334422414635 
0895598897 470774007 
-0856224485318131325 
-0808855721934550922 
-0753934869239738285 
-0691974694940161476 
-0623553678524653054 
-0549310594426269679 
-04699385046 10241705 
-03861781477 18139676 
-0298810459167464775 
-02086460901 76033601 
.0116484483922677346 
.0018939393939393939 


-07758792403848083 16 
-0773546 125373055485 
-076656081 1933174706 
-07549653 10459179728 
-0738829357473896064 
-0718249996197639705 
-0693350992870948951 
-0664282092328427060 
-0631218117275497954 
.0594357916635436557 
-0553923169189143105 
-0510157049479865135 
-0463322763486516145 


-0251176983936010865 
-0193532985662516799 
-01347207 10944706066 
-0075066587628515885 
-0012195121951219512 


-064785433 1007237352 
-0646495667764058446 
.0642425376733883324 
06356605301 11866414 
.0626229501980326849 
-0614171849294469270 
-0599538145961051888 
-0582389770704910239 
-05627986496 1 1358469 
-0540846954421352817 
-0516626757839646535 
-0490239647292254072 
-0461796298735572735 
-0431416012270164702 
039922621 1435904434 
.036536 1908133766348 
.0329965 135057250782 
-0293184347098430027 
-025517379167 1367047 
.0216092842408220077 
.017610526816783 1063 
.0135378293815818172 
-0094080469628555945 
-0052366002997777452 
.0008503401360544218 
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TaBLeE 1—Continued 











Xkn akn 
n = 65 
-0000000000000000000 -0487112532912551590 
-0486919954825551174 .0486534844338217974 
-0972684989276217033 .0484803 148828198298 
- 1456142922319649050 .0481921553771416351 
.1936147045111101818 .0477896893990385782 
.2411558840858965205 .04727387 15529625068 
- 2881250685259474757 -0466459253013279525 
-3344108521095384325 -0459073400625635379 
.3799034500654639962 .0450598676783326357 
.4244949589701384082 .0441055182582980979 
-4680796126822756626 .0430465554 122800447 
-51055403320807 20304 .041885490881 1053505 
-5518174759018253800 .0406250785788684 140 
-5917720684203420792 . 0392683080607 104407 
.6303230428642674229 .0378183974315733463 
-6673789605555873413 .0362787857 126874349 
. 7028519289179370194 .0346531246837989863 
- 7366578099449523957 -0329452702203 182977 
.7687164197616275394 -0311592731456411222 
. 7989517 188043673806 .0292993696198098245 
.8272919921669404206 .02736997 10863165681 
-8536700196814493087 .0253756537989475783 
-8780232353249457 100 .0233211479495077525 
- 90029387556 16070109 .0212113264134911180 
-9204291162429907785 -0190511931200282741 
.9383811976826655570 .01684587 10220649483 
.9541075374600222128 -0146005895504222135 
- 9675708302669 136143 -0123206711194442981 
.9787391331919988895 .0100115149738883044 
-9875859307695091355 .0076785701335481211 
.9940901501184231212 .0053272417215879641 
- 998235858985 1681587 .00296203254 12562160 





1.0000000000000000:900 
, (4) 
é+1) _ (4) P,.s(th (i) 
(4) Ukn ae pp’? 7.(4)\ ein A. 
I n—1(Lkn 


Since [9] 


(5) (1 — 2°)Pa(x) = n[Pra(x) — xP,(2)] 
and 
(6) (1 - t’)Py (x) = 2rP’,(x) — n(n + 1)P,(2), 


A reduced to the following formula which only involves Legendre polynomials: 
ies - Pra) — 22 Pralal? 
Qakn [Pro(tin) — tin Paa(tin )] — nPra(ain 

The iteration was terminated when |A| became smaller than 2-“. The initial 
values xi) were taken to be the arithmetic means of two successive zeros of Fai) 
tabulated by Davis and Rabinowitz. [10, 11] With this choice of initial value, 
convergence was achieved in at most six iterations. 

The values in the table were checked on the computer by calculating the sums 
rm Gintin for r = 0, 2, 4, 8, 16, 32. For n such that r < 2n — 3, these sums 
were equal to (2/r + 1) to 19 decimal places. In addition > 7 __m azn was computed 
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by hand from the final manuscript and was found to deviate from the correct 
value 2 by at most 4 units in the 19th decimal place. This accuracy is acceptable 
since we were summing rounded numbers. A further check was made on the com- 
puter by calculating >°7-. zi. and []%;2:.. Each sum checked to 19 decimal 
places with the true value (n — 2)(n — 3)/2(2n — 3) while the square of each 
product checked to at least ten significant figures with the true value 


n+1 
("ti) m+n 


These two identities are derived as follows: 


(8) P,3(2) = ges Z. (—1)° n—1\f2n — 2 — 2t Pe eel 
t=0 t n—1 
‘ _ 2 — « a 2 
P,-(2) = 7 1) > (—1)* (" 7 * 2 ; ‘ — Pp 2t) Pata 
t=0 _ 


= _ In —2—2 
=o) «>, (-(" ‘ _ fe N(x — 1 —2t)(z*)”*" 


(9) t=0 3 1 
x [ay (2°)”* + a,(2*)”* +--+ + apo + anil 


MT II (x + Xen) (X — Len) = Ag x II (x* — aj, ). 
k=1 = 


n—1\/2n —4 
- as, 1 Oe gta 


ao n—1\f/2n —2 Ta 2(2n — 3) 
( 0 Y( n— ac 1) 
and 


oe 1 eee et 
Bi an (9)? St og _m+i1\m+1 


Hence 


M 3 
&, 








> 
ll 
= 





k=l A n—1\f/2n — 2\, 2 fan—-2° 
(chee Ey 
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A Method for the Numerical Evaluation of Finite 
Integrals of Oscillatory Functions 


By I. M. Longman 


1. Introduction. In two previous publications [1, 2] the author has demon- 
strated a method, based on Euler’s transformation of slowly convergent alternat- 
ing series, for the numerical evaluation of infinite integrals of oscillatory func- 
tions; this can be used in many cases by a double application for the evaluation 
of finite integrals of oscillatory functions. For example the integral 


(1) [ Hoax 


where f(z) may have a very large number of oscillations in the range of integration 
can conveniently be evaluated as 


(2) [ se a - [ 4) dz. 


However in physical problems the finiteness of the range of integration is often 
associated with a kind of natural boundary of f(z), such that it is impossible to 
extend f(x) to values of x beyond the upper limit b while preserving the general 
character of f(z). Analytically speaking, s = 6 may be a branch point of f(z). 
Alternatively, it may be possible to extend the range of integration to infinity as 
in equation (2), but the infinite integrals may not converge. As an example of the 
branch point difficulty we can consider the integral 


(3) [= [ (a’ — 2*)' sin x dz, 
0 


which, if a is large, would be very difficult to compute by straightforward numerical 
integration owing to the large number of oscillations, and we clearly cannot apply 
the method of equation (2). An instance of this kind of difficulty in a physical 
problem is given in Pekeris [3] where the solution of a problem of sound propaga- 
tion in a layered liquid is given in the form of infinite integrals where the character 
of the oscillatory integrand changes abruptly at a certain point in the interval of 
integration. 

The present paper gives an extension of Euler’s transformation and applies 
it to the numerical computation of integrals such as (3). 


2. Description of the Method. For simplicity let us start by considering alter- 


nating series, and suppose that we wish to calculate the sum of a series of the 
form 


(4) S = 1 — 1% +2 — --- + (—1)"0, 
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where the »,’s are all positive and slowly decrease in numerical value as 7 increases 
from 0 to n. For example such a series might be 

(5) S = 9999' — 9998’ + 9997) — --- + 17. 


In order to transform (4) into a form convenient for computation let us con- 
sider the associated power series, 


(6) S(x) = v9 — vy + vex” — --- + (—1)"0,2", 
which reduces to (4) when x = 1. If we multiply by x and add we obtain 
(1 + x)S(x) = 1 — (% — v)x 
+ (vo — v,)a° — +++ + (—1)"(vn — Ona)” + (—1)"0,2"", 
or, with the usual notation for differences, 
Av; = Vign — V5 


+1 
A’ —— A’vi4n - A’0; ’ 





we have 
(1 + x)S(x) = 10 — (Avp)a + (Av, )2” — +>» + (—1)"(Avg_)a” + (—1)"0,2"7?. 
This gives 
' vo + (—1)"v,2"" 
S(z) = 
(7) " itz 


— y[Avo — (Ar)x + (Av2)a* — --- + (—1)""(Av,4)2”"], 


where y = x/(1 + x). A second application of this transformation to the series 
in square brackets in (7) yields 


vo + (—1)"vn2™"* Avo + (—1)" "(Avy)" 





= i+z i+a " 

+ y'[A’ro — (A’r:)x + (A’v2)a® — +++ + (—1)"*(A°en-2) 2" ], 
and p applications give 
S(x) = (vo — ydvo + y’ Arvo — +>» + (—1)? "y” A”) /(1 + 2) 


+ (—1)"[vnv"™ + (Avgps)2"y + (A’vn2)2""y? + ++: 
+ (A vn-pyi)a” Py? "(1 + 2) 
+ (—1)?y"[A?v0 — (A?r,)a + (A?r)2” — --- 
+ (-1)” ?(A7v,_»)x” ”], psn. 
Putting x = 1 we have as a transformed form of (4) 
S = [(1/2)vo — (1/4) Avo + (1/8) A’ — --- + (—1)?"2°7(A? "w)} 
(8) + (=-1)"[(1/2)on + (1/4) Avna + (1/8) A70n-2 + +++ + 277A? dnp 4] 
+ 2°-°(—1)?[A?u — A”, + Avg — +--+ + (—1)” 7A? vn_5], psn. 
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This result is of course exact, but for large values of p (assuming that the high 
order differences are small) the later terms in the first two square brackets and 
the whole of the third square bracket in equation (8) can be neglected since 2” 
will be negligible. 

Assuming then, that n is large, we have for (4) the excellent approximation 


S = (1/2)vo — (1/4) Aro + (1/8) A’ — --- 
+ (—1)"[(1/2)on + (1/4) Av, + (1/8) A’vn-2 + --- | 


which represents a kind of double application of Euler’s transformation. 


3. Examples. The utility of the series (9) will be demonstrated by a number of 


examples. 
1. Consider the series 
(10) S=1-1/24+ 1/3 — 1/44+--- — 1/1000. 


To evaluate (10) by means of equation (9) we split off the first eight terms (whose 
differences do not decrease very rapidly) and evaluate 


(11) S’ = 1/9 — 1/10 + 1/11 — --- — 1/1000. 
The contribution at the first eight terms of S is 0.634524, and so 
S = 0.634524 + S’. 


The differencing of the first few and the last few terms of S’ is shown in Table 1. 
From this, by means of equation (9), we obtain 


S’ = (1/2) K 0.111111 + (1/4) XK 0.011111 + (1/8) X 0.002020 
+ (1/16) X 0.000505 + (1/32) X 0.000156 + (1/64) X 0.000057 
— (1/2) X 0.001000 + (1/4) X* 0.000001 = 0.058123. 


Thus S = 0.692647. As a check we can calculate (10) as the difference between 
two infinite series 


S=(1-1/2+1/3-i/4+---) 
(12) — (1/1001 — 1/1002 + 1/1003 — --- ) 
= In2 — (1/1001 — 1/1002 + 1/1003 — --- ). 


The series in parentheses in (12) can be evaluated by applying Euler’s transforma- 
tion, and its sum is easily shown to be 


0.000499. 
Thus working with equation (12) we obtain 
S = 0.693147 — 0.000499 = 0.692648, 


agreeing with our previous result. 
2. In example 1 we were able to extend our series to infinity so that we really 
had no need to use the transformation (9). We now consider, however, an example 





TABLE 1 




















a at as 
97 = 0.111111 
—11111 
10 = 0.100000 2020 
—9091 — 505 
11+ = 0.090909 1515 156 
—7576 — 349 —57 
12+ = 0.083333 1166 99 24 
— 6410 — 250 —33 
13+ = 0.076923 916 66 
— 5454 — 184 
14> = 0.071429 732 
— 4762 
157 = 0.066667 
994-1 = 0.001006 
—1 
995 = 0.001005 
—1 
996 = 0.001004 
—1 
997—' = 0.001003 
—1 
998-1 = 0.001002 
—1 
999 = 0.001001 
—1 
1000 = 0.001000 
TABLE 2 
a? ; at * 5 x 
9999? = 99.995000 
— 5000 
9998! = 99.990000 —2 
— 5002 +2 
9997! = 99.984998 0 1 
— 5002 +3 
9996? = 99.979996 3 
— 4999 —6 
9995! — 99.974997 —3 
— 5002 
9994! = 99.969995 
15? = 3.872983 
— 131326 
144 = 3.741657 —4780 
— 136106 — 563 
13? = 3.605551 — 5343 — 122 
— 141449 — 685 —35 
12? = 3.464102 — 6028 —157 —27 
— 147477 — 842 —62 
11? = 3.316625 — 6870 —219 
— 154347 — 1061 
10! = 3.162278 —7931 
— 162279 


9! = 3.000000 
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of a series which cannot be so extended. Such a series is given in equation (5). 
This series has 9,999 terms and its direct evaluation would be a rather lengthy 
computation. However it is easily computed by an application of our transforma- 
tion in equation (9). Differencing the first few terms, and a few terms near the 
end of the series we obtain Table 2. We split off the last eight terms of the series 
(which do not difference so well) and obtain 


S = (1/2) X 99.995000 + (1/4) x 0.005000 — (1/8) x 0.000002 
+ (1/2) X 3.000000 — (1/4) x 0.162278 — (1/8) x 0.007931 
— (1/16) X 0.001061 — (1/32) x 0.000219 — (1/64) x 0.000062 
—-#47-6¢+45 -— 443! — 24 1! = 50378853. 


3. Now let us turn to the evaluation of oscillatory integrals. We will illustrate 
this by evaluating the integral (3) for a = 100x7. We have 


1007 
(3’) [ = [ (100°x” — 2x”)! sin x dz. 
0 
By splitting up the range of integration, J can be expressed as the series 
99 = 
(13) > (-1)* | (100°? — (rx + 2)"} sin x dr 
r=0 0 


which is of the form (4). A few integrals near the beginning and near the end of 
the series (13) were evaluated on an IBM 709 computer using a 16 point Gaussian 


integration formula. The results are tabulated and differenced in Table 3. Applying 
our method we obtain 


I = (1/2) X 628.30915 + (1/4) 0.06285 — (1/8) X 0.06284 
+ (1/16) X 0.00004 — [(1/2) X 325.85292 — (1/4) & 10.14092 
— (1/8) X 0.41146 — (1/16) X 0.03382 — (1/32) & 0.00472 
— (1/64) X 0.00089 — (1/128) X 0.00022] + 315.26077 
— 304.17027 + 292.52472 — --- — 60.96022 

= 298.43558. 


As a check we apply the method differently to Table 3, obtaining 
I = 628.30915 — 628.24630 + (1/2) X 628.12061 — (1/4) X 0.18857 
— (1/8) X 0.06298 — (1/16) X 0.00002 + (1/2) X 238.71325 
— (1/4) X 14.76162 — (1/8) X 0.96185 — (1/16) x 0.14230 
— (1/32) X 0.03311 — (1/64) X 0.00997 — (1/128) X 0.00362 
— (1/256) X 0.00150 — 222.79836 + 209.46181 — --- — 60.96002 
= 298.43557 


agreeing with our previous result. 











a ug 


1 628.30915 — 
2 628.24630 
3  628.12061 
4  627.93204 
5 627.68049 


6 627.36594 
79 381.13726 
80 372.75841 
81 364.07840 
82 355.07512 
83 345.72330 
84 335.99384 
85 325.85292 
85 325.85292 
86 315.26077 
87 304.17027 
88 292.52472 
89 280.25486 
90 267.27464 


91 253.47487 


92 238.71325 
93 222.79836 
94 205.46181 
95 186.30583 
96 164.30583 
97 139.47917 
98 108.12528 
99 60 . 96022 


— 6285 
— 12569 
— 18857 
— 25155 


— 31455 


— 837885 
— 868001 
— 900328 
— 935182 
— 972946 
— 1014092 


— 1059215 
— 1109050 
— 1164555 
— 1226986 
— 1298022 
— 1379977 


— 1476162 
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TABLE 3 
a? At As 
— 6284 
—4 
— 6288 
—10 
— 6298 
—2 
— 6300 
— 30116 
—2211 
— 32327 —316 
— 2527 —67 
— 34854 — 383 —22 
— 2910 —89 
— 37764 —472 
—3382 
—41146 
— 49835 
— 5670 . 
— 55505 — 1256 
— 6926 —423 
— 62431 — 1679 212 
— 8605 —635 
— 71036 — 2314 362 
— 10919 —997 
—81955 —3311 
— 14230 
— 96185 


—150 
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4. Conclusion. An extension of Euler’s transformation has been presented and 
its use in the computation of finite integrals of oscillatory functions demonstrated. 
It is believed that its application will render feasible the numerical solution of 
physical problems hitherto regarded as intractable. 
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Products of Laguerre Polynomials 


By Joseph Gillis and George Weiss 


1. Problems occasionally arise in theoretical physics where one wishes to ex- 
press the product of two linear combinations of Laguerre polynomials as a linear 
series of these same polynomials. The purpose of this note is to investigate the 
coefficients C,,, where the following definitions are used: 


(1) L,(x) = > (—3y (") x’ /r! 


r+s 


(2) L(x) L(x) = > Crat L(x) 


t=|r—s| 


The limits of the sum in (2) follow, in fact, from repeated application of the re- 
currence formula for the L,’s. This can be written [1] in the form 


(3) aL,(x) = —(n + 1)Dngi(x) + (2n + 1)L,(2) — nLg1(z). 


It follows from the orthogonality properties of Laguerre polynomials that 
(4) Cru = [ e*Le(2) Lax) Lia) de 
0 


and, in particular, is symmetric in r, s, t. A closed formula has been obtained for 
Cy by Watson [4]. We begin by obtaining the same formula by a very simple 
argument. In §3 we derive a simple recurrence relation suitable for rapidly gen- 
erating the coefficients as needed when working with a high speed computing ma- 
chine. This will be seen to be more useful in practice than the formal expression in 
equation (7). 

2. It is known [2] that the Laplace transform of L,(x) is p ““(p — 1)" while 
that of L,(x)L,(2) is 


h+k\ (p — 1)" .. P(p — 2) 
( h ) 2a a - ee ° 


Hence, taking Laplace transforms of both sides of (2), we get 


r+es\ +41 rts . —- . MP — 3) 


= > Crp (p — 1)*. 





With gq = 1 — p’ we have 


(6) > Cag = (" . ‘) q*Fil—r, —s; —r — 8; (2g — 1)/q'l. 
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Comparing coefficients of q‘ in (6) gives us 


Cus site > (r + i n)! (—3)”"""* (, n ) 


» ni(r — n)i(s — n)! mn—-r—s+t 





(7) bp )! 
om i utW ns r s—n)! 
-i-o Xu 2 (r — n)\(s — n)\(2n — p)p — n)! 
where we have written p = r + s — t. This is equivalent to Watson’s formula. 
The limits of the sum in (7) are defined by the requirement that none of the 
arguments of the factorials can be negative. It is easily confirmed from this, in- 
cidentally, that there can be no terms in the sum if ¢ lies outside the range (|r — s |, 
r+s). 
3. To establish a recurrence relation we base ourselves on the well-known result 
[3] that, if u(z), v(z) satisfy the normalized differential equations 


u”(x) + I(x)u(xz) = 0 
v(x) + J(x)o(z) = 0, 
then y = wv satisfies the equation 


d fy" +2704+J)y¥ + ('4+Jd"%)y) = 
(9) al TJ |+U-Dy =0 





(8) 





provided that J # J. In case ] = J, y satisfies the third order equation 
(10) yy” + Aly’ + 2I’'y = 0. 


Applying this result to the differential equation satisfied by Laguerre polynomials, 
after normalization, we obtain the following equation for y = L,(z)L,(x)(r 2 s) 


D(y) = 2’y°” + 2(5 — 4x)y’” + [4 + (20 — 15)x + 52'ly” 


il 
on + [(30 — 8)—4(o — 3)a — 22°ly’ + [8 — 3(¢ — 1)4+2(o — 1)zly = 0 
where 
o=r+s8s+1 
(12) : =r-— 8s. 


In fact, equation (11) holds whether or not r = s. We now substitute from (2) 
into zD(v), making use of the following formulas: 


(13) aly = —(1—2)L/ — th, 
(14) 2h,” = (2— (¢+2)¢4+ ob! + (2 -—2)ke, 
(15) aL? = —{6 — 2(2t+ 3)x + (2t + 3)a° — ahh? 

—46—(¢+3)2+ 24k. 


Of these equations, (13) is simply the differential equation of LZ, , while (14) and 
(15) are obtained from it by differentiation followed by substitution from (13) 
itself. We thus obtain, after some reduction, 


(16) 2D(L.) = r2(1 — 22)L/ + {Qra* + afe — ¢ — r(2t + 3)RLr 
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where 

(17) r=o-—-t-—-l=r-+s-t 
But, [1], 

(18) aL! = (¢+ IL — (t+1—2)Ky. 


Substituting from (18) into (16) and making repeated use of (3) leads finally to 

aD(Li) = 2r(t + 1)(t + 2)Lege — (t + 1)[r(4t + 5) + (8 — P)LDias 
+ (2+ 1)[(¢+ 1)r + % — OL, — (8 — #) Lia. 

For fixed r, s write C,,, = Az. It follows from (19) that 

(20) xD( 2) Ad,) = >> Bi, 


where 


(19) 


B, = —(¢+ 1)[8 — (¢+ 1)74en 
(21) + (2¢+ 1) — + (t+ 1)o — (t+ 1)"JA, 
— ifs — (t — 1)? + (44+ 1)(6 — tJ Aci + 2t(t — 1) (0 —t + 1)Ace. 
We thus obtain the recurrence relation 
(t+ 1)[8 — (t+ 1)'Aen = (264+ DP -— P+ (E+ 1)(o —t — DIA 
(22) — 8 — (t— 1)? + (44+ 1)(6 — t)JAra 
+ 2(t —1)(o —t+ 1)Ace. 


We know that A; = 0 fort < 5. However (22) becomes indeterminate for ¢ = 6 — 1 
and so A; has to be calculated independently. This was to have been expected since 
equation (11) is homogeneous. It follows immediately from (7) that 


(23) A; = Ceaen — (") 


When working with an electronic computer it will nearly always be more efficient 
to use (22) and (23) than (7). All that one need store is the binomial coefficients 
(23) and a sub-routine for effecting (22). 

4. It may be of interest to consider some values of C,,, for special values of 
r, s, t. The value of C,,,,--. is given by (23) and we deduce, by means of (22), that 


Cr.e1-0141 = — 28 (") and 
(24) ’ 
bee 3s 10 + DP = 80 (’). 


We can calculate directly from (7) that 


Caiasae = q ; *) 


(25) © .0,040-8 = 2(r + o 1) (" + en 4 


r—l 


Crisrpo-2 = (278s —r —s +1) (" +s— ‘). 


r—l 





f 
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One might also draw attention to some elementary arithmetical properties of 
the coefficients C,,;. In the first place, it is clear from (7) that the sign of C,,. is 
that of (—1)” and so of (—1)"****. Again, it follows from the same formula that 
C,.t is always an integer. For the expression 


(r+s—n)! 
(r — n)(s — n)! (2n — p)! (p — n)! 
is an integer, being a multinomial coefficient. Also, the occurrence of the term 


(2n — p)! in the denominator imposes the limitation 2n = p on the range of n. 
The result is immediate. Finally, we deduce, by putting xz = 0 in (1) and (2) that 


(26) > Crt = 1 





and, from symmetry, that the summation in (26) could equally be taken over s or r. 

5. The relation (22) is, as we have said, well adapted to machine work. It is 
rather complicated for hand computation and the authors are indebted to a referee 
who drew their attention to the alternative relation 


(27) (r + 1)Cr41,2.¢ = (t + 1)C,.0,t41 + 2(t aa 7)C+.2.t sad Wicains + i yab-2 . 


This follows immediately from the orthogonality and recurrence relations of the 
Laguerre polynomials, and is very much simpler arithmetically than (22). There 
is no doubt that many other relations of this type could be found. However, for 
machine computation they would all share the disadvantage of (27), namely, the 
increased programming complications involved in varying two of the subscripts. 
They would also be slower to generate since, to arrive at a given r, s, t, one would 
have to progress through a much larger number of intermediate terms. 
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TECHNICAL NOTES AND SHORT PAPERS 


The Formula-Controlled Logical 
Computer “Stanislaus” 


By F. L. Bauer 


The evaluation of a formula of propositional calculus is considerably simplified 
if this formula is written in the parenthesis-free notation of the Warsaw School, 
[1]. The Warsaw notation may be formulated in the following way: 

There are symbols for operations, e.g.: N for negation; K for conjunction; 
A for disjunction; E for equivalence; C for implication; and symbols p, q, r, s, ¢ 
for variables. 

A variable is a formula. A formula preceded by the symbol N is a formula. 
Two juxtaposed formulas preceded by any one of the symbols, K, A, E, C are a 
formula. Evaluation of such a formula is done in the following way: Each of the 
variable symbols p, g, r,--- has a value 0 or 1. The operation symbol acts on the 
value of the one or two formulas governed by it giving the value of the compound 
formula. 

It was remarked in 1950 by H. Angstl, [2], that a mechanical evaluation of a 
formula, written in Warsaw notation without brackets, can be done in the fol- 
lowing easy way: Each of the variable symbols is represented by a box with one 
output, the negation by a box with one input and one output, and the other opera- 
tion symbols by a box with two inputs and one output. The meaning of a formula 
in the Warsaw notation is given by Angstl’s rule: The first input of each opera- 
tion symbol is to be connected with the output of the next following symbol, 
either variable or operation. The symbol N excepted, the second input of each 
operation symbol is to be connected with the first remaining free output of a 
symbol going from left to right. This may be demonstrated by an example, which 
uses Stanislaus present capacity of eleven symbols: the tautology of transitivity of 
the implication [(p — q) & (q—1r)] — (p — 1), written in Warsaw notation 


CKCpqCaqrCpr 


and represented according to Angstl’s rule by the diagram Fig. 1. 

Following the Angstl rule, there is an easy way to build a mechanism for the 
evaluation of a logical formula, especially for a switching network which does the 
necessary connections of the logical units corresponding to the boxes mentioned 
above. Of course, these units can be built up in the usual way, e.g., with relays. 
But the proper switching may be done automatically by punching the formula on 
a keyboard. 

We found building such an apparatus highly instructive on account of our 
interest in formula-translation, because it is a model of a computer which imme- 
diately obeys formulas in a common language. This means that, apparently in 
contrast to the apparatus of Kalin and Burkhart, [3], and later modifications, logical 
formulas are directly “‘written” (e.g., typed on a keyboard) rather than wired 
manually by connecting outputs and inputs on a switchboard. 

It seems quite clear that this idea of direct formula control is not restricted to 
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output eS input 


second input 


HE G82 Ga AGB 


Fic. 1.—Example of a mechanical evaluation of a formula. 











logical operations and the values 0, 1 only. Indeed, our interest is primarily devoted 
to arithmetical formulas, and we have obtained (in collaboration with Dr. K. 
Samelson, of Mainz) some results we believe to be remarkable. They will be pub- 
lished in the near future. 

The wiring design of Stanislaus was done in December 1950. With the friendly 
help of colleagues, the toy was built in the following years, mostly from surplus 
material. The work suffered long interruptions, since it had a very low priority in 
our working program. Stanislaus was finished at the end of 1956 and presented in 
a lecture by H. Angstl, “‘Vorfiithrung eines logistischen Rechengerites fiir den 
Aussagenkalkiil,” held on January 8, 1957, in the Logistisches Seminar der Uni- 
versitat Miinchen. Stanislaus is now sometimes used for demonstration purposes. 
There might be possibilities of serious use outside the scope of our interest. 

In the language of compriter engineering we constructed, influenced by the ma- 
terial provided to us, something like a parallel-in-formula computer, but serial-in- 
fortaula operation would have been possible too, and the essential idea would have 
been the same. Indeed, the Burrough Truth Function Evaluator as described in 
1954 by Burks, Warren and Wright, [4], is realized in this way. 

For reasons of simplicity of the design we were going as far as providing each 
column in the keyboard with its own logical unit. The logical units are built up 
conventionally by relays. Also the truth values are represented by zero or working 
voltage on the output of the five switches corresponding to the values 0 or 1 of 
the variables p, q, v, s, t. Our interest was not devoted to these points, but to an 
automatic connection of the logical units by a mechanized switching system which 
works according to Angstl’s rule. This switching system contains an additional 
feature which shows whether the string of symbols on the keyboard is a formula or 
is not well-formed. 

Details may be seen from the wiring diagram, Fig. 2. Application of the formula- 
controlled logical computer is most simple. One only has to push buttons on the 
keyboard which operates like a desk calculator keyboard and therewith to “write 
in’ the formula. The truth value of the variables involved may be set on switches, 
A red or yellow light flashes on according to the value of the formula being “wrong” 
or “right’’ for these values. A special switch provides the test whether the string 
of symbols ‘written in” is a formula (blue light) or not well-formed. 

The keyboard scheme may be seen from Fig. 3. The formula written-in is 
KNpNq or pA@ in common notation. The switch at the bottom of the key- 
board is in position “Check for meaningfulness.” After changing it to the other 
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Evaiuotion —-® Verification 


Fic. 2.—Wiring diagram. 

Pushbuttons on the formula keyboard, columns 0-10; O, blank; N, negation; C, conjunc- 
tion; D, disjunction; EZ, equivalence; J, implication; p; qg, corresponding; r, variables; s; ¢. 

Relays, columns 0-10; X, Y, logical circuits; x, y, corresponding contacts; A, B, trans- 
mission lines; a, b, corresponding contacts. 

Relays and switches, left side: p, q, r, s, t, switches for truth values input. 

Evaluation-Verification switch for evaluation or check; K, evaluation relay; k, corre- 
sponding contact. 
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Fig. 3.—Keyboard scheme. 
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side, the yellow result light will flash on, indicating the truth value of KNpNq 
for p = wrong and g = wrong according to the left position of the switches for 
the variables p and gq on the right side of the keyboard. 
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Evaluation at Half Periods of Weierstrass’ Elliptic 
Function with Rectangular Primitive 
Period-Parallelogram 


By Chih-Bing Ling 


The purpose of this paper is to evaluate the following Weierstrass’ elliptic func- 
tion at half periods [1], 


(1) &: = Y(w), €2 = M(w), €: = Yes), 
where 2; and 2w: are double periods of the function and a; is defined by 
(2) w + ow + w = 0. 


This paper tabulates only the values of the function whose primitive period- 
parallelogram is a rectangle with 2w, = 1 and 2w2. = ai, where a = 1. 
The three functions in (1) form a set of distinct roots of the cubic [1] 


(3) x — pr —q=0, 
where 
(4) p = ldo, g = 3506, 
and 
=! 1 
oH = 





m,h=—20 (2ma; + 2nwe)* 


x 1 x x 1 
=22) 5+22 2 


m2 i mo (m + nai)” 


(5) 





The accent on the summation sign denotes the omission of simultaneous zero values 
of m and n from the double summation. 
The cubic (3) indicates that 
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(6) €: + @ + e; = 0. 


Also, since ¢; , ¢2 and es are distinct, the discriminant (4p° — 27q°) of the cubic 

does not vanish. As will be seen later, in the present case both o, and o¢ are real 

and (4p° — 27q°) is positive. This implies that all the roots of the cubic are real. 
The evaluation of o, and gg is facilitated by using the known relation [2] 


(7) cot = ++ > ( i - +) 


o \mr + x mr 





where the accent on the summation sign denotes the omission of the zero value 
of m from the summation. By repeated differentiation of Equation (7) and sub- 
stitution of ix for x, it is found that 

















= 1 a 2 4 1 
(8) m=—o(mr-+ix)* 3sinh’r  sinh‘z 
< 1 ah 2 to +. od 
mano (me +ix)®  15sinh¢z sinh‘z — sinh*z’ 
Hence we have 
Loa Sa 
* 45 " 3 sinh?xa 
(9) aati 
= doe > ae 
945 15sinh?za 
where 
K, = sinh*ve > ( a. ) 
(10) nai \sinh’nra 2 sinh‘nza 








K, sinh?ra 5 ( ee Sap yeee ) 


“=i \sinh’nza = 2 sinh‘nrza_~—= -2 sinh®nza/ © 


Consequently, we find 
(11) 4p’ —27q¢ 5Ki+7K2 , 100K, — 147K, , 2000K,’ 


1672 3 sinh?7ra sinh‘7ra sinh*ra © 





With the aid of known tables [3, 4], values of K,, Ke, and then o,, os and 
(4p° — 27q°)' are computed to 16D for a = 1(0.25)2(1)6 and © as shown in 
Table 1. 

The subsequent evaluation of ¢; , ¢2, and e3; requires the solution of the cubic 
(3). It appears that one of the roots, e; , can be easily evaluated to 16D as shown 
in Table 2 by using Newton’s method or otherwise, but difficulty arises in eval- 
uating the other two roots for in most cases they are almost equal. However, they 
can be separated by forming a new cubic 


(12) a+px—q =0 


whose roots are the differences of the roots of the cubic (3). Let (e, — es), (e2 — é3) 
and (es; — e,) be the roots of the new cubic. We have 


, 


(13) = (€; — &2)(€2 — €3) + (€2 — es) (es — 1) + (€3 — &1) (1 — 2) = —S3p, 
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(e; — €2)°(e2 — e3)*(es — &)” = 4p — 27q°. 
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Consequently, by taking a positive sign for q’, the new cubic is in the form 
(14) a — 3px — (4p — 27q°)' = 0. 
From this cubic, values of (e: — e;) and then e: and e; are computed to 16D as 
shown in Table 2. 
It is mentioned that the values of the function, for 0 < a < 1 or in general for 


w./w, purely imaginary, can be computed from the tabulated values with the aid 
of the following relation [1] 


(15) @(Az | Nor, Awe) = Nz | wr , we) 


where J is a constant, real or complex. 

The writer wishes to express his thanks to Mr. C. P. Tsai for his assistance in 
performing the numerical computations. The writer also is deeply grateful to 
Professor C. W. Nelson of the University of California, Berkeley, for checking the 
manuscript and verifying all the numerical values in Tables 1 and 2 by independent 
calculations. Thanks are also due to the referee of the paper, who suggests a different 
method of computation [5] without solving the cubic equation. 
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A Note on the Nonexistence of Certain Pro- 
jective Planes of Order Nine 


By Raymond B. Killgrove 


1. Introduction. Every finite projective plane may be coordinatized in at least 
one way |1]. In this process some line is chosen to be the line at infinity, and the 
points not on this line are represented by an ordered pair of elements. The elements 
x and y for any point (2, y) on a given line of the plane satisfy the equa- 
tion y = x-mob, where m and b are specific elements for the given line. This ternary 
operation on x, m, and b includes an additive loop in a special case. 

A sequence of SWAC computer routines has been written to search for all 
planes having a specific additive loop in an appropriate ternary ring. Using these 
routines, a complete search had been made previously using the elementary Abelian 
group for the additive loop [2]. Now a complete search has been made using the 
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cyclic group for the additive loop. No planes were found. This entire note parallels 
the paper [2]. 


2. The lines consistent with the cyclic additive pencil. It is sufficient to find the 
affine planes of order nine where the lines parallel to the line y = x are determined 
by the cyclic group. In particular these lines are of the form y = x + c where y 
may be computed for each x and a given c by addition modulo nine. 

We represent any possible line for a geometry as a permutation on nine marks 
in the following way: a number a and its image b under the permutation define a 
point (a, b) of the line. All of these lines can interesect in at most one point with the 
lines of the cyclic groups. We can also restrict ourselves to the pencil of lines going 
through the origin (0, 0). Only 225 lines were found which satisfy these criteria. 

By applying the permutations defined by the additive loop, the full list of 
2,025 lines consistent with the cyclic group may be found. Having the lines x = ec, 
y = c, in addition to y = x + c, one needs to find sets of 63 lines from the 2,025 
which satisfy the affine plane axioms. 

To save time one considers the automorphisms which preserve the 27 known 
lines, and in particular the subgroup of automorphisms which fix (0, 0). The cyclic 
group of order nine has such an automorphic subgroup of order 6. When this sub- 
group is applied to the 31 lines of the 225 containing the point (1, 2), one obtains 
25 conjugate classes of lines. The conjugate classes are placed in some preference 
order. In this case the order chosen was 4, 16, 8, 9, 10, 11, 13, 15, 18, 19, 20, 23, 
24, 25, 2, 3, 6, 7, 17, 21, 12, 14, 1, 22, 5. All but the last five classes mentioned in 
this order have six members. Conjugate classes 12 and 14 have three members, 
1 and 22 have two members, and 5 has only one member. 


3. Construction of pencils through (0, 0). Using a card sorter one finds the lines 
of the reduced set consistent with each line through (1, 2). Conjugate classes of 
lines preceding the current sorting line in preference order are removed also. This 
set of sorted cards is given the computer, which in turn does further sorting and 
the printing out of complete Latin squares. On the average each conjugate class 
yielded eight Latin squares. These Latin squares represent possible pencils of lines 
through (0, 0). 

Next we use a set of 225 transformed lines which are obtained by applying one 
of the non-identity permutations to the original 225 lines. The computer then lists 
the transformed lines consistent with a given Latin square. In the previous ele- 
mentary Abelian group case [2], whenever the list of transformed lines associated 
with some Latin square could themselves form a second Latin square, one usually 
finds a geometry. In this cyclic group case a second consistent Latin square was 
never found. Hence the search ended here, as there are no possible planes in this 
situation. 
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A Note on Rational Approximation 


By Robert W. Floyd 


It is suggested by plausible reasoning and confirmed by experience that the 
error of an nth degree polynomial approximation, in the Chebyshev sense of least 
maximum error, to an analytic function, is roughly a multiple of the n + Ist 
Chebyshev polynomial, 7',,:(z), on the interval of approximation. Therefore if 
the nth degree polynomial f*(z) is equal to the function, f(z), on the roots of 
Tn41(x), we expect that f*(x) will be a satisfactory approach to a Chebyshev 
approximation of f(x). 

Because f(x) is analytic, it may be represented with negligible error in the 
interval of approximation by a polynomial p(x) of sufficiently high degree; e.g., 
a truncated Taylor’s or Maclaurin’s series. Applying the division algorithm for 
polynomials, 


p(x) Qo(2)*Tr4i(x) + ro(x) 
Trai(x) = q(x)-ro(r) + (2) 
ro(a) = go(x)-ri(x) + r2(x) 


r(x) qs(x)-re(x) + 713(x), ete., 


where the degrees of the r; form a strictly decreasing sequence. From these equations 
we may write r;(x) = a;(x)-p(x) + b:(x)-Tnii(x), where a; and 6; are defined 
recursively by 


Q; = Qi-2 — GiGi, a_, = 0, a s=1 
b; = bie — git bia ’ bi = 1, b_s = 0. 


It may be proven that the sum of the degrees of a;(x) and r;(x) is at most n. The 
first set of equations may be written p(x) = [r;(x)/a:(x)] — [bi(x)/a:(z)]-Trai(z), 
so that r;(x)/a,(a) is a rational approximation to p(x), exact wherever 7',4:(2x) 
vanishes. Since 7',4:(2) < 1 in the interval of approximation, b;(x)/a,;(x) provides 
a bound for the error of the approximation. If b;(2)/a,(a) is nearly constant on the 
interval of approximation, the error oscillates between n + 2 extrema of nearly 
equal magnitude, and the method of approximation is justified, for Chebyshev 
approximation is characterized by an error which oscillates at least n + 1 times 
between positive and negative extrema of equal magnitude. For the particular 
case 7 = 0, a; = 1, and 7(2) is a polynomial approximation to f(x) of degree at 
most 7. 


For example; f(z) = & = 1+ a+ (2°/2!) + (2°/3!) +--:; 
p(x) =1+2+ 52° + .1666 66672° + .0416 6667z* + .0083 33332" 
+ .0013 8889x° + .0001 9841x' + .0000 2480z° + .0000 02762". 


For —1 < x S$ 1,| p(x) —f(x) | $ 3.0 X 10°". Tr(a) = 642° — 1122° + 562° — 72. 
Then go = (317.5625 + 38.752 + 4.31252") x 10°; 
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ro = 1 + 1.0000 22232 + .5000 02712" + .1664 89132° + .04164497z" 


+ .008686592° +- .0014 32292°; 
| p(x) — ro| = | qo|-| T(z) | $ 3.61 XK 10° (—-1 S21). 
Therefore | f(z) — ro| < 3.91 X 10°° (—1 S x S 1). Dividing T;(z) by 70, 
q: = —270,998.81 + 44,683.688z. 
r, = 270,998.81 + 226,314.152 + 90,815,4582° + 22,832.3912° 
+ 3,846.3890z* + 381.20482°. 


do = 1; bo = — 
is —q 5b = 1 + Qo 
Therefore 
p(x) — oo 41 Fae p(,), 
ay ay n n 


The second term on the right is 


.1394 0940 + .0368 865982 + .0056 2810542" + .0019 269840z° 
— 270,998.81 + 44,683.688z 


whose absolute value is bounded by 8.121 X 10°’ for —1 < x < 1. Thus & may 
be approximated on this interval by 


1 + 8351 11232 + .3351 13862" + .0842 52742° 
‘. - + .0141 93382" + .0014 06672" 





T(x) 


n «1 = .1648 85182 
where the error is bounded by +(3 X 10°’ + 8.1 X 10°’) = 41.1 X 10°° 


Armour Research Foundation 
Illinois Institute of Technology 
Chicago 16, Illinois 


1. C. Hastines, Approximations for Digital Computers, Princeton, 1955, p. 47-64. 
2. F. B. HrtpesRanp, Introduction to Numerical Analysis, McGraw- mil Book Co., New 
York, 1956, p. 389-395. 


3. NBS Seseamp Matuematics Serres, 9, Tables of Chebyshev Polynomials S,(z) and 
C,(x), U. S. Govt. Printing Office, Washington, D. C., 1952, p. 16-18. 


The Complete Factorization of 2'” + 1 


By K. R. Isemanger 


The integer 2 + 1 is divisible by 2“ + 1 = 17-353-2931542417 and the 


quotient, 2% — 2“ + 1, is divisible by 241-7393. There remains the formidable 


problem of factoring the resultant quotient N, where N is the integer 
1 73700 82040 22350 83057. 


Received July 3, 1959; in revised form, September 17, 1959. 
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By a method of exclusion, using small prime moduli, the following binary repre- 
sentations of N were obtained: 


88073 81316" + 98046 34351” 
82086 25547" + 2°-3°-11°-37"-4243”-7-17-19-73 


= 
I 


= 1 28399 72408’ + 7°-11°-13°-3137°-3-19-79-199 
= 23334 17999° + 2*-13°-17°-2707°-3-7-79-89-199 
1 06383 10009" + 2”-17°-3853°-3- 107-167-257 
= 1 05791 40907? + 2°-3°-13°-3169"-17-23-29-89- 167 


= 1 58709 07595" — 2°-3°-17°-37°-3061°-7-13-29 


2=zaaz2a2a22aaa 
Il 


= 1 47403 24637? — 2°-7°-23”-737-2803"-3-7-239 
51299 86183" — 2°-7°-: 


K 
= 
Il 
oe 
| 
J 
— 
~ 
—— 
bo 
ou 
© 
b 
ae 
wo 
“J 
~ 
ae 
— 
—) 
~J 
i) 
ad 
Jo) 


= 
I 


49963 46199° — 2°-3*-11°-17°-2711°-13-23-257 
If these ten relations are written as congruences in the form 
« = Dy'(mod N) 
and then are multiplied together, there results the congruence 


A* = B’(mod N) 
where 
A = 3030 76720 24193 56872 


B 


85638 82032 43137 98848. 


The greatest common divisor of A + B and N was found to be 9 86182 73953, 
which yields the factorization 


N = 17613 45169-9 86182 73953. 
The primality of both factors was verified on the SILLIAC computer at the 
University of Sydney. 
This result supplements information previously published by R. M. Robinson [1]. 


19 Snape Street 
Kingsford, Sydney 
New South Wales, Australia 


1. RapHakEt M. Rosrnson, “Some factorizations of numbers of the form 2" + 1,’’ MTAC, 
v. 11, 1957, p. 265-268. 
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REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 


1{B, C, K, S].—A. H. Wapsrra, G. J. NuGu & R. van Liesuout, Nuclear Spectros- 


copy Tables, North Holland Publishing Company, Amsterdam, 1959, viii + 
136 p., 27 em. Price $8.90. 


This book consists of a well-selected collection of tables and graphs of interest 


to nuclear physicists. A slightly condensed list of contents follows: 


Chapter I. Mathematical Data. 
Common logarithms. Powers of 10 and 2, for exponents 1.00(0.01)9.99, 
4D. Cube roots of integers 1(1)999, 5D. Brief discussion of least-squares 
method, and a graph for checking the consistency of least-squares com- 
putations. Table of Gaussian distribution and its integral for arguments 
0.00(0.02)3.48, 4D. 

Chapter II. Tables of Atomic Constants. 

Chapter III. Elements and Isotopes. 
Names, atomic numbers, and symbols of elements. Atomic weights and 
abundances of isotopes. 

Chapter IV. Heavy Particles. 
Range energy curves for protons, deuterons, and alpha particles. Straggling 
of protons in air (graph). Magnetic rigidity of protons, deuterons, and 
alpha particles (table). Half-lives for alpha disintegration and spontaneous 
fission of heavy elements (graphs). 

Chapter V. Electrons. 
Range energy curves for beta-particles. Saturation backscattering coefficient 
(graph). Magnetic rigidity of electrons (table). Discussion of shapes of 
continuous beta-spectra. Reduced Fermi function f(z, p) and beta-decay 
functions L and M for negatrons and positrons (tables). Screening correc- 
tion for negatrons and positrons (graphs). Discussion of beta-decay transi- 
tion probabilities and their values. Ratios of electron capture in different 
electron shells. Computation of log ft values (nomogram and graphs). 
K capture/positron emission ratios (table and graphs). 

Chapter VI. Gamma Rays. 
Half-thickness of some substances for absorption of gamma rays vs energy 
(graph). Photon absorption cross-sections vs energy (table). Energy of 
Compton scattered gamma rays (discussion and graph). Gamma-decay 
half-lives (discussion, nomogram, and graphs). 

Chapter VII. X-rays and Auger Electrons. 
Electron binding energies in different shells for all elements (table). Rela- 
tive intensities of K X-ray components and of K-Auger electrons (table). 

Chapter VIII. Angular Distributions and Correlations. 
Brief discussion of angular distributions and correlations involving photons, 
alpha particles, beta rays and K conversion electrons. Tables of F,(L, L’, 
ji, J) coefficients. 

Chapter IX. Nuclear Models. 
Discussion of nuclear mass formula, nuclear shell model, collective model, 
magnetic moments, quadrupole moments, and gamma- and beta-decay 
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probabilities. Graphs of Nilsson level scheme of single particle orbits in 
spheroidal potential. Table of measured ground state spins of odd-A and 
odd-odd nuclei. Tables of Clebsch-Gordan coefficients. 

Chapter X. Calibration Standards. 

Tables of standard gamma and electron lines and of standard alpha rays. 
Gamma-ray absorption coefficient in Nal crystals. Table of standard nu- 
clides for calibration of gamma-ray spectrometer. 

The tables and graphs have been presented so as to be easily read, and the 
quality of the printing is good. Much of the material is used frequently by nuclear 
physicists but is widely scattered in the literature. Thus, this book should prove 
very helpful to people in the field of nuclear physics, and this reviewer recommends 
it highly. 


MicHet A. MELKANOFF 
University of California 
Los Angeles, California 


2[D, L].—L. K. Freven, J. W. Turney & D. R. Perersen, Seven-Place Table of 
Iterated Sine, The Dow Chemical Company, Midland, Michigan, 1959. De- 
posited in UMT File. 


Following a detailed description of the method of computation employed, the 
authors give a 7D table of the nth iterated sine function of x for n = 0(.05)10, 
and x = k(2/20), where k = 1(1)10. It is stated that the computations were per- 
performed on a Datatron 204, and the results are considered correct to within 
5-10". 


J. W. W. 
3[G]._K. M. Howe tt, Revised Tables of 6j-Symbols, U. Southampton Math. 
Dept., Research Report 59-1, 1959, xvi + 181 p., 33 cm. 


The Wigner 6j-symbol has been defined by Wigner in general in connection with 
the reduction of the triple Kronecker product of any simply reducible group. In 
these tables this group is taken to be either the three-dimensional rotation group 
or the two-dimensional unitary group. The symbols are denoted by 


Sis ja i 
\ hike ks 

where the quantities 7, , --- , ks are integers or half-integers. If we let 

Jo=fAitptis, A=Atht+k, Jz=jtht+ks 

Js=jotht+kh 
Ki =~jtjpe +h t+ ke Ke = jit jst hi + ks 
K; = jotjsthe+hks, 

then the explicit expression for the 6j-symbol is 


(RAR = Uw - soul. + on 


. > (—1)'(¢ + 1)!/{I] (K, — ): I] @ -— J.)!}. 
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These symbols are invariant under the 144 products of separate permutations 
of the K, alone and the separate permutations of the J, alone. The symmetries of 
the 6j-symbols are made use of in the organization of the tables. In the present 
(revised) edition of the tables a wider group of symmetries is exploited than in the 
earlier edition. 

In the tables the 6j-symbols are classified in terms of a set of six ordered param- 
eters, and the tables are arranged in descending “‘speedometer”’ order of these param- 
eters. Rules for determining the values of these parameters from given j; , --- , ks 
are included. 

The square of the value of the 6j-symbol is printed in the table, together with 
its correct sign. In addition the entries are written as rational fractions in terms of 
their prime factors. Thus the fraction —4/5+/21 is written as —2-2-2-2/3-5-5-7. 
A composite member, all of whose factors are greater than 103 is printed as if it 
were a prime. 

The tables were duplicated from stencils cut directly from the output tape 
of a Ferranti Pegasus Computer. The program used in the computer was essentially 
the same as that used for the earlier version of the tables. It is stated that, “It 
seems reasonably sure that these tables are as reliable as, and more comprehensive 
than, the first set of tables of the 6j-symbol.” 


A. H. T. 


4{H, Mj.—H. J. Gawuix, Zeros of Legendre Polynomials of orders 2-64 and weight 
coefficients of Gauss quadrature formulae, A. R. D. E. (Armament Research and 
Development Establishment) Memorandum (B) 77/58, Fort Halstead, Kent, 
December 1958, 25 p. 


1 
The Gauss quadrature formula is [ f(x) dx = Dr, As(a,) with A,, a, 
1 


chosen so that the equality is valid whenever f(z) is any polynomial of degree 2n — 1 
or less. The quantities a, are the zeros of the Legendre polynomial P,(z), while 


1 * P, (x) 
P,/(a,) Liz — a, 

The memorandum here reviewed gives 20-decimal values of a, and A, for all 
zeros a, of each polynomial P,,(a), for n = 2(1)64, and is by far the most extensive 
and accurate of such tables available. 

The most extensive of the previously published tables is that of Lowan, Davids, 
and Levinson [1], which gives a similar table to 15 decimals for n = 2(1)16. This 
table has been reproduced several times, in whole or in part [2]. It has been supple- 
mented by a table by Davis and Rabinowitz [3], which gives 20-decimal values for 
n = 2, 4, 8, 16, 20, 24, 32, 40, 48. Discrepancies between Davis & Rabinowitz and 
Gawlik amount to only a unit in the 20th decimal in several places, with the proba- 
bility that Gawlik is the more accurate; both tables would thus appear to be reliable. 


J.C. P. M1tier 


A, = 





The University Mathematical Laboratory 
Cambridge, England 


1. Amer. Matn. Soc., Bull., v. 48, 1942, p. 739-743; v. 49, ae. B 939 (errata). 

2. See, for example, NBS Applied Mathematics Series, No. 37, Tables of Functions and of 
Zeros of Functions, 1954. 

3. NBS, Jn. of Research, v. 56, 1956, p. 35-37. 
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5{1].—-A. N. Lowan, The Operator Approach to Problems of Stability and Conver- 
gence of Solutions of Difference Equations and the Convergence of Various Iteration 
Procedures, Scripta Mathematica, New York, 1957, x + 104 p., Office of Tech- 
nical Services, Department of Commerce, Washington 25, D. C. 26 cm. Price 
$3.00. 


In the numerical solution of linear second-order differential equations by differ- 
ence methods, one has to solve 


(1) Aur a Buu _ Cyte-r = D,. 


where A; , B,C; , Dy are sparse matrices with regular structure, the u’s are vectors, 
and the integer subscripts refer to time-steps, iteration cycles, etc. In many impor- 
tant cases A; , B, , C, are independent of k and closely related. Stability and mesh- 
convergence of “stepping-ahead” solutions (parabolic, hyperbolic equations) and 
convergence of iterative solutions (elliptic equations) can be discussed in terms of 
the operators A; ‘B, and A, ‘C;. In the von Neumann technique, as formalized 
and extended by the reviewer, one essentially takes the Fourier transform (in the 
extended sense) of (1), thus introducing immediately the eigenvectors and eigen- 
values of these operators. This is equivalent to a change of coordinates in vector 
space under which A; , B; , C; take very simple forms. Various authors, however, 
have retained the original coordinates and worked directly with (1); prominent 
among these are 8S. Frankel, D. Rutherford, A. Mitchell, P. Lax, R. Richtmyer, J. 
Douglas, J. Todd, and (unpublished) C. Leith. Professor Lowan here gives a de- 
tailed, connected account of this second method which he calls the “operator ap- 
proach”. He covers the usual parabolic, elliptic, and hyperbolic partial differential 
equations, homogeneous and non-homogeneous, with some attention to equations 
with variable coefficients. He discusses stability and mesh-eonvergence for parabolic 
and hyperbolic equations, and convergence of the various iterative schemes for 
elliptic equations. Original contributions, besides the organization of the material, 
include a discussion of iteration schemes for solving ‘‘implicit” difference approxi- 
mants to parabolic and hyperbolic equations and a novel ‘“‘second-order”’ Richard- 
son scheme for elliptic difference equations. In addition, Professor Lowan has writ- 
ten down a number of “folk theorems’’, rather widely known but unpublished. 
There are eight “Sections” and six ‘‘Appendices”. Most of the typographical errors 
have been detected by the author and listed on the ‘Errata’ sheet. However, the 
figures on page 5 and at the top of page 33 should be corrected; in the second line 
from the bottom of page 44, read (—1)"*'4/2 sin rhx/M; page 26, line 5 and page 
47, line 18 should show the respective ouclifications ‘““r > 1” and ‘‘r < 3”. More- 
over, the reader should be aware that the estimates of truncation error in difference 
solutions given during the discussion of mesh-convergence require that various func- 
tions be “sufficiently continuous.’’ Those of us interested in the numerical solution 
of partial differential equations are indebted to Professor Lowan for a very worth- 
while addition to the literature of this field. 
Morton A. HyMAN 


IBM Research Laboratory 
Yorktown Heights, New York 
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6{K].—Cuar.es W. Dunnett, “Tables of the bivariate normal distribution with 
correlation 1/+/2,” 1958, 28 em. Deposited in UMT File. 


The bivariate normal probability distribution, 


re 1 ry f - (z — 2rey + v) 
L(h, k; r) ev gay & exp| 3 ‘Tr. dz dy, 


has been tabulated by Karl Pearson [1] for r = —1.0(.05)1.0. The present tables 
were prepared to avoid the necessity of interpolating in Pearson’s tables when 
r = +1/+/2. This case arises in certain double sampling procedures in which 
probability statements concerning X and X + Y jointly are required, where X and 
Y are independent normal chance variables with the same variance. 

The tables were computed on a Royal McBee LGP-30 electronic computer, by 
numerical quadrature, using the relation 


Lh, k;r) = [ E - r(: - 2) | (2) dx 


where f(x) = (1/+/2z) exp (—2°/2) and F(x) = / f(x) dx. The function is tabu- 


lated for r = 1/+/2 in Table I and for r = —1/+/2 in Table II for positive values 
of its arguments, h varying in steps of 0.1, and k varying in steps of 0.1+/2. All 
entries are given to six decimals and should be correct to this number of places. 

The function can be determined for negative values of its arguments by using 
the relationships 








L(—h, k;r) = L(—~#,k;r) — L(h, k; —r) 
L(h, —k;r) = L(h, —~;r) — L(h, k; —r) 
L(—h, —k;r) = L(h, k;r) + 1 — L(h, —«&;r) — L(-—~,k;r) 


where L(h, —©;r) = L(—«,h;r) = 1 — F(h), which is the right-hand tail area 
of the univariate normal distribution. In order to avoid the necessity of consulting 
tables of F(h), these values are included in the table. 

Tables III, IV and V were computed from Tables I and II using these relation- 
ships. The error in the entries in Tables III and IV should be no greater than a unit 
in the sixth decimal place. The error in the entries in Table V should be no greater 
than two units in the sixth decimal place. All tables have been deposited in the Un- 
published Mathematical Tables repository. 

AvuTHoR’s ABSTRACT 


1. Karu Pearson, Tables for Statisticians and Biometricians, Part Il, Cambridge Uni- 
versity Press, 1931. 


7(L]._E. A. Cusrova, Tablitsy funktsit Besselia ot detstvitel’nogo argumenta i 
integralov ot nikh (Tables of Bessel functions of real argument and of integrals 
involving them), Izdatel’stvo Akademii Nauk SSSR, Moscow, 1958, 524 p. + 
loose card, 28 em. Price 45 rubles. 


This important volume in the now familiar series of Mathematical Tables of the 
Computational Center of the Academy of Sciences was initiated by V. A. Ditkin. 
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The main table occupies pages 23-522 and gives values for 
xz = 0(.001)15(.01)100 
of the eight functions 
Jn(z), Jina) = f ‘Yaa 





¥,(2), Yig(a) = f 2 an, 

where n = 0 and 1. The values are to 7D or, near singularities at the origin, 7S 
Auxiliary functions (detailed below) are provided for x = 0(.001).150. Over the 
range x = 1.350(.001)15 there are no differences, and linear interpolation provides 
results correct within two units of the last place. Second differences are required for 
some functions in parts of the ranges x = .150(.001)1.350 and x = 15(.01)100; 
the quantities printed (when greater than about 16) are sums of two consecutive 
second differences, for use in Bessel’s formula. 

The values of Jo(x) and J;(x) are rounded to 7D from the well-known Harvard 
tables, and are included for convenience; the values of the other six quantities 
result from original computations, and may be checked only partially against 
previous, less extensive, tables which are listed in a bibliography. The integrals 
Jig(x) and Ji;(x) were computed by Simpson’s rule on an electronic computer and 
other machines. The functions Yo(x), Yi(x), Yio(x), and Y%,(x) were evaluated 
on the electronic computer BESM, using Taylor series and asymptotic expansions. 
All values were checked by differencing. By means of formulas given on pages 11-12, 
the integrals of Jo(u), Ji(u), Yo(u), and Y;(u), not divided by u, may be simply 
expressed in terms of the tabulated functions. 

The nine auxiliary functions given on pages 17-19 are all tabulated for x = 
0.(.001).150 to 7D without differences. The functions are: 


Tig(x) = Jio(x) + In $x 
Co(z) = Yo(x) — (2/r)Jo(x) nz E(x) = (2/r){ Jto(x) + In 2} 
Ci(z) = a{¥i(x) — (2/e)Ji(x) naj} E(x) = (2/e){Jnx(x) — 1} 
Dox) = (2/2) Jo(x) "o(x) = Yio(x) — (2/x) In a{ Jin(z) + 3 In 2} 
Dy(x) = (2/r)Ji(2) F\(x) = 2{ Ya(x) + (2/r) In afl — Ji(z)} 


It is stated that the errors do not exceed 0.6 final unit, except that they may 
attain one final unit in the case of the auxiliary functions Co(x), Ci(x), Fo(x), and 
F,(2) ° 

A table of the Bessel coefficient }¢(1 — ¢) for ¢ = 0(.001)1 to 5D without differ- 
ences is given on page 523 and also on a loose card. 

A. F. 


8[L].—Orro EMERSLEBEN, “‘Werte einer Zetafunktion zweiter Ordnung mit Argu- 
ment s = 2,” Bearbeitet in der Abteilung Angewandte Mathematik der Uni- 
versitat Greifswald, Greifswald, 1956. 
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This paper contains a tabulation of the function 
C) eritetis) 


u(2,y) = > > = 


k=—oo l=—oo RP+E + P 


where the prime denotes the fact that the term with k = / = 0 is omitted. 

The table is given for values z and y in the range 

s2>x2y20, 

where x and y = 0(0.01)0.5. The entries are given to six and sometimes seven deci- 
mal places, and are said to be accurate to at least two units in the last decimal place. 

In the calculation of this table, use was made of the seven-place tables of the 
exponential integrals published in 1954 by the U.S.S.R. Academy of Science, Insti- 
tute of Exact Mechanics and Computation. 

A. H. T. 


O[L].—HerMAN H. Lowe x, Tables of the Bessel-Kelvin Functions Ber, Bei, Ker, 
Kei, and their Derivatives for the Argument Range 0(0.01)107 .50, Technical Re- 
port R-32, National Aeronautics and Space Administration, Washington, D. C.., 
1959, 300 p., 26cm. 2. 


These tables provide an elaborate and attractively arranged compilation of 
decimal! values of the Bessel-Kelvin functions (frequently referred to simply as the 
Kelvin functions) of the first and second kinds of order zero, together with their 
first derivatives. Approximations to ber and bei and their first derivatives appear in 
floating-point form to generally 13 or 14 significant figures. On the other hand, the 
number of significant figures given for ker and kei and their first derivatives vary 
from 9 to 13, according to a pattern explained in the detailed introduction, which 
also describes the construction of these tables and the checks applied to the tabular 
entries. The calculations were performed on an IBM 650 calculator using the Bell 
Telephone Laboratories Double-Precision (16-figure) Interpretive System. 

In addition to the checks applied by the author, the reviewer collated the values 
of ber x and bei z with similar data given by Aldis [1] to 21D, for the range x = 
0.1(0.1)6.0. No discrepancies were detected. 

The range, precision, and accuracy of the tables under review establish them as 
the definitive tables of the Kelvin functions at the present time. 

J.W. W. 

1. W. SreapMAN A pis, ‘‘On the numerical computation of the functions Go(z), G:(z), and 

J(rVi),” Roy. Soc. London, Proc., v. 66, 1900, p. 32-43. 


10(L].—NumericaL CompuTaTion Bureau, Report No. 11, Tables of Whittaker 
Functions (Wave Functions in Coulomb Field) Part 2, The Tsuneta Yano Me- 
morial Society, 1-9 Yuraku-cho, Chiyoda-Ku, Tokyo, Japan, 1959, 52 p., 26 
cm. Price $3.00. 


The first part of these tables was reviewed in MTAC, v. 12, 1958, p. 86-88. The 
earlier review contains some errors and fails to give complete information. 
The functions considered are defined thus: 


Gi. + iF: = exp (-3 :) exp | -i G r— «)| W it.t4sy2( —tz). 
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The earlier review has 3 in place of / + 3 in the subscript of the Whittaker function, 
W, and incorrectly uses G;:,1/2(2) and F¢,1;2(2) for the case 1 = 0. Both the previous 
review and the present tables fail to define o, = arg '(/ + 1 + iz). 

The table now reviewed is concerned with values of F:,0(2) for large E. Values 
of this function are presented to five decimal places, together with 6,’ and 6,’, where 
k= - and r = 2x/(2k), corresponding to r = 0(.1)10 and k = 0(.01)1, that is, 


er 


ae is observed that for k = 0, that is, > ~, 


Fro(z) = V2i(2V/2r); Geo(z) = V2r¥,(2V/2F), 


which are identifiable as Bessel-Clifford functions multiplied by 2r. This relation 
allows comparison with appropriate data in a publication [1] of the National Bureau 
of Standards; such a comparison has revealed 24 errors (all due to rounding) in the 
tables under review, thereby suggesting a general accuracy therein within one unit 
of the fifth decimal. 


J.C. P. Mrtuer 
The University Mathematical Laboratory 
Cambridge, England 


1. NBS Applied Mathematics Series, No. 28, Tables of Bessel-Clifford Functions of Orders 
Zero and One, U. 8. Government Printing Office, Washington, D. C., 1953. 


11[L, M].—W. W. Gersss, G. E. Reynotps, M. R. Hoss, & C. J. Drang, Jr., 
Table of S(x) and its First Eleven Derivatives, Vol. 1, 2, 3, Air Force Cambridge 
Research Center, Bedford, Massachusetts, 1958, 27 cm. 


The tabulated function S(x) defined by 


2 ° 
. & 

z{ sin 3 

S(x) _ [ u/2 du 


is related to the sine integral Si(x) by 
S(x) =2 | Sica) “ 1— 003, 


For ease in computation in the design of antennas, the function S(xz) and its 
first eleven derivatives are tabulated to six decimal places for 


x = 0°(1°)18,000°. 


The introduction gives the characteristics of the functions, reduction formulas, 
power series representations, asymptotic expressions, integral representations, 
differential equations, transforms, addition formulas, etc., and the method of 
computation. 

The tables were computed using the IBM 650 calculator. 


IRENE A. STEGUN 
National Bureau of Standards 


Washington 25, District of Columbia 
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12([L, M].—R. Hensman & D. P. Jenkins, Tables of? & [ e dt for Complez z, 
T z 


Royal Radar Establishment, Malvern, Worchestershire, England. Deposited 
in UMT file. 


¢ 


The function 2 / edt has been tabulated to 6 decimal places for 


0(0.02)2.00 in the real part of z and 0(0.02)4.00 in its imaginary part, and also for 
0(0.1)10.0 in both real and imaginary parts. Second differences are given and sec- 
ond-difference interpolation in the appropriate table gives 6-decimal accuracy for 
the whole range covered. The error in using linear interpolation need not exceed a 
unit in the fourth decimal. 


AutTHors’ SUMMARY 


13{L, M].—A. V. Hersney, Computing Programs for the Complex Exponential 
Integral, NAVORD Report No. 5909, NPG Report No. 1646, U. S. Naval 
Proving Ground, Dablgren, Virginia, June 1959, iii + 16 p. Figures and tables. 
27 cm. Astia Document Service Center, Armed Services Technical Information 
Agency, Arlington Hall Station, Arlington 12, Va. 


In this report there appears a detailed discussion of the use of asymptotic series 
with remainder in the evaluation of the complex exponential integral on the Naval 
Ordnance Research Calculator (NORC). Brief descriptions are also given of two 
NORC subroutines for the calculation of the exponential integral and the sine and 
cosine integrals of a real argument. 

A series expansion of the rernainder of the asymptotic series is presented, and is 
used to evaluate the remainder with improved accuracy. 

The author also describes the construction of a rational approximation to the 
exponential integral that is valid over the negative half of the complex plane outside 
the unit circle. In appended tables appear approximations to 13S of the coefficients 
of two polynomials of the fourteenth degree whose quotient gives values of the func- 
tion ze “Ei(z) to within a maximum relative error in the absolute value of 2-2 
10”. 

The appendix also includes a table of values to 13S of the real and imaginary 
parts of Hi(z), corresponding to x = —20(1)20 and y = 0(1)20. These results 
were obtained on the NORC by means of double-precision arithmetic, using sixteen- 
point Gauss integration over each unit interval, beginning with z = — 100, where 
the value of the integral was considered to be negligibly small. 

Comparison by the reviewer of these data with corresponding entries in an ex- 
tensive earlier set of tables [1], carried to 6D and 10D, revealed only three in- 
stances of rounding errors in the latter, all of such size as to lie within the guaran- 
teed limit of a unit in the last decimal place. 


i - & 
1. NBS pogiet Mathematics Series, No. 51, Tables of the Exponential Integral for Complex 


Arguments, U.S. Government Printing Office, Washington, D. C., 1958. See also MTAC, v. 13, 
1959, p. 57-58. 
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14{L, M].—K. A. Karpov, Tablitsy funktsit F(z) = | é dx v kompleksnoi oblasti 
0 


(Tables of the function F(z) = [ é* dx in the compiex domain), Izdatel’stvo 
0 
Akademii Nauk SSSR, Moscow, 1958, 518 p. + 2 inserts, 27 cm. Price 61 rubles. 


This is a companion volume to the tables reviewed in MTAC, vol. 12, p. 304- 
305, and completes the tabulation of the error function in the complex plane. The 


present volume contains 5D or 5S values of the real and imaginary parts of the 
function 


F(z) = [ dr = u + w 
4 


for z = pe”,0 < p < mo, x/4 < 6 S x/2 and 6 = O. The quantity po depends on 
6 and decreases from po = 5 for @ = 2/4 to po = 3 for 6 = 2/2. An exception is 
6 = 0, for which po = 10. In the introduction, a diagram is given representing the 
intervals in 6 and the value of po for each @, and a table indicates the intervals in p 
in various parts of the volume. As in the earlier volume, the diagram is reproduced 
on a cardboard inset, which serves also as an index to the numerical tables. 

The introduction gives integral representations and series expansions for u and 
v, graphs of u and v as functions of p for selected values of 6, relief diagrams of u 
and v over the sector of tabulation, a description of the tables and numerical ex- 
amples showing their use, some useful numerical values, values of cos 26, sin 26, 
and values of (2n + 1)@in radians for n = 0(1)5 and for those values of @ included 
in the tables. There is also a one-page auxiliary table of t(1 — ¢)/4for0 <¢ < 1. 
This table, together with a nomogram for finding A’t(1 — t)/4, where A? designates 
the sum of two consecutive second differences for use in Bessel’s interpolation 
formula, is reproduced also on a cardboard inset. r 

Using the symmetry properties of F(z), this function can now be evaluated on 
the real axis and in a sector of half-angle 45° to both sides of the imaginary axis. 
Between them, Karpov’s two volumes contain a very satisfactory tabulation of 
the error function in the complex plane. 


A. ERDELYI 
California Institute of Technology 


Pasadena, California 


15[L, M, S].—Y. Nomura & 8S. Katsura, “Diffraction of electric waves by circular 
plate and circular hole,” Sci. Rep. Ritu, B-(Elect. Comm.) 10, No. 1, 1958, 43 p. 


The problem of the diffraction of a plane electromagnetic wave by an infinitely 
thin, perfectly conducting, circular disk of radius a, and the problem of the diffrac- 
tion of such a wave by a circular hole of radius a in a plane conducting screen are 
discussed. The method of solution involves the expansion of the two-component 
Hertz vector in terms of hypergeometric polynomials. The solution is valid for all 
frequencies. However, convergence is pooi when ka = 27a/d becomes large. Tables 
of values are included of the real and imaginary parts of 





Gr, = (2m + 4 + 1) [ Sas ae dé 
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for m = 0(1)4, n, v = 0(1)5, and ka = 0(.25)5, where J,(£) is the Bessel function 
of the first kind of order p. The function X%, is defined in terms of the equations 


2 C82 @4,. 

v=0 
Tables for X%, for n = 0(1)4, v = 0(1)2, 5, m = 0(1)4, over the same range of £, 
are also included. 

In addition, tables are given which are useful for the calculation of the field 
distribution at large distances, and tables are given which enable one to determine 
the current distribution on the plate and the electric field distribution on the whole. 

All table entries are given to four decimal places. However, no indication is 
given as to the numerical method of evaluating the table entries nor as to their 
actual accuracy. 

A. &P. 


16{P|.—Cuartes J. THorne, Temperature Tables: Part 1. One-Layer Plate, One- 
Space Variable, Linear, NAVORD Report 5562, U. S. Naval Ordnance Test 
Station, California, 1957, iv + 711 p., 28 em. 
This table is concerned with listing the solution of the heat conduction equation 
in a plate of finite thickness, with heat transfer at both faces. In mathematical form 


Ux = kU, 0< xX < L, t> 0 
KUx = —h(U; — U) X = 0, t>0o0 
KUx = —h(U — Uo) X = L, t>0 
U = Up ¢t = 0, 0O<X<L 
where the conductivity K, density p, specific heat c, diffusivity h = cp/K, heat 
transfer coefficients h; and ho , and stagnation temperatures U; and U» are assumed 
to be constant. L is the plate thickness, z is the distance, and ¢ is the time. For the 
tables a dimensionless system of variables is adopted, ie., x = X/L, kL’r = t, 


u = (U — U»)(U; — Uo), ao = hol /K, a; = h,L/K. Then the above problem 
becomes 


Mao = MU; S< 2 <.1, 7>0O 
Uz = —a,(1 — u) z = 0, tT>0O 
Uz = —aou z= 1, 7>0O 
u=0 r= 0, o<2< 1 


The analytical solution of this problem is 


ao(1 -+ ar) bes —B,21 
2 an, 
a + ai + aa; >> » D’(B,) 


-{a,[8, sin B, — a» cos B,] sin 8,2 + a {8, cos 8B, + a sin 8,] cos 8, 7}, 


ulz,r) = 1— 





where 


D’(B) = —B(2 + a + a;) sin 8B + (—B + ao + a; + aoa;) cos B, 
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and the 8,’s are the positive roots of 
(aoa; — 6°) sin B + B(ao + a;) cos B = 0. 


The table gives the dimensionless temperature u, where 0 < u < 1, for the follow- 

ing ranges of parameters 

x = 0(.01).02(.03).05(.05).3(.1).7(.05).95(.03).98(.01)1 

.001 (.0005 ) .002(.001 ).008( .002 ).01(.01 ).08(.02).1(.1).8(.2)1(1)8(2)10(10) 

80( 20) 100( 100 )800( 200 ) 1000 

a; = .001, .002, .004, .006, .01, .02, .04, .06, .1, .2, .4, .6, 1, 2, 4, 6, 10, 20, 40, 60, 
100, 200, 400, 600, 1000. 

a = 0, .001, .004, .01, .04, .1, .4, 1, 4, 10, 40, 100, 400, 1000. 

The tabular entries are given to five figures, with better than three of the 
figures being accurate. For the most part, the error appears to be one or two units 
in the fifth figure. 

The table should be very useful to those people who are engaged in design work 
involving heat transfer, as, for example, rocket nozzle design. The introduction 
also contains a generalization of the heat conduction problem defined above which 
can be solved by means of the tables. 

The tabular entries are printed with reasonable clarity. There are, however, a 
few obvious misprints in the introduction. 


ll 


+ 


R. C. Roserts 
Naval Ordnance Laboratory 


Silver Spring, Maryland 


17(S|.—-G. A. BartHotomew «& L. A. Hiaes, Compilation of Thermal Neutron 
Capture Gamma Rays, Report AECL No. 669, 1958, 146 p., 27 em. Available 
from Scientific Document Distribution Office, Atomic Energy of Canada 
Limited, Chalk River, Ontario, Canada. Price $2.50. 


This report presents a compilation of energy, absolute intensities, and spectral 
distribution of gamma rays produced by capture of thermal neutrons, together with 
a complete bibliography of information on this subject through June 1, 1958. The 
results obtained from measurements at the Chalk River Laboratories over several 
years using the pair spectrometer have been reviewed. These results have been modi- 
fied where necessary such that all intensity determinations are presented on a uni- 
form basis. The accuracy of pair spectrometer intensity measurements is discussed. 

Included in the tables are the energies and intensities (photons per hundred 
captures) of resolved gamma rays obtained from experiments in which absolute 
intensities were determined. References to other data not tabulated is also given. 
Where an appreciable portion of the gamma ray spectrum is unresolved, a spectral 
distribution curve is given. Most of the curves plot the number of gamma rays per 
capture per Mev as a function of energy. Results published by the Moscow group 
are included. A rough measure of accuracy and completeness of the results is given 
with each tabulation and with most of the curves. 

K. SHURE 
Westinghouse Electric Corporation 


Bettis Atomic Power Division 
Pittsburgh, Pennsylvania 
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18[S].—M. E. Rose, Internal Conversion Coefficients, Interscience Publishers, Inc., 
New York, 1958, xxii + 173 p., 31 cm. Price $6.25. 


When a nucleus is in an excited state for which the excitation energy is insuffi- 
cient for the emission of nuclear particles, the de-excitation will proceed predom- 
inantly by either one of two competing mechanisms. Either a y-ray will be emitted 
or the nuclear excitation will be transferred to one of the orbital electrons resulting 
in the ejection from the atom of this electron. The latter process is referred to as 
internal conversion. If N, is the number of conversion electrons emitted per second 
and N, is the number of photons emitted per second, the internal conversion coeffi- 
cient @ is defined as 


a = N./Nyg. 


This book gives a comprehensive account of a ten-year program devoted to the 
calculation of internal conversion coefficients. It contains a ten-page introduction 
which gives a precise and thorough account of the physical and numerical approxi- 
mations made in the course of the calculations of the tables. These in turn consti- 
tute the bulk of the work, 164 pages. 

The conversion coefficients are strongly dependent on k, where kmc’ is the tran- 
sition energy, Z the atomic number, L the angular momentum change, and on Az, 
the parity change. The tables list values of a, and 6, , (L = 1, 2, 3, 4, 5), the co- 
efficients for 2” electric pole and 2” magnetic pole conversions, respectively, fork = 
0.05(0.05 )0.2(0.2)1.0(0.5)2.0 and Z = 25(1)95. 

Also included in the tables are values of certain radial matrix elements R;(m) 
and R,(e) for the K shell. These are uncorrected for screening and for finite nuclear 
size. 

The author lists the sources of error in determination of the radial wave func- 
tions, which are a fundamental set of intermediate quantities in the calculation of 
the tables, and expresses the view that all of these errors are small and amount at 
most to 1-3 per cent. He expresses the view that the irreducible minimum error 
involved in any calculation of internal conversion coefficients, aside from nuclear 
structure effects, is just smaller than the experimental error in the best measure- 
ments now available. 

The author states that, “interpolation in the tables will be necessary only in the 
energy variable k. For this purpose interpolation on a log a or log 8 versus log k is 
advisable since the plots of the conversion coefficients on a log-log scale show very 
little curvature.” 

a. a Be 


19[S].—C. H. Wesrcort, Effective Cross Section Values for Well-Moderated Thermal 
Reactor Spectra, Report AECL No. 670, 1958, 27 p., 27 em. Available from Scien- 
tific Document Distribution Office, Atomic Energy of Canada Limited, Chalk 
River, Ontario, Canada. Price $1.00. 


This report is concerned with the determination of an effective neutron absorp- 
tion cross section, which cross section is recommended for calculations involving 
reaction rates. The effective cross section is defined in terms of a neutron density 
distribution per unit velocity. The neutron spectrum assumed consists of a Max- 
wellian distribution at a temperature T°K plus an admixture of a 1/Z distribution 
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of flux per unit energy interval, the admixture being controlled by an epithermal 
index r. For r = 0, the spectrum is a pure Maxwellian. 
Using this spectrum, the effective cross section, ¢ can be written as 


o = o220(g + rs) 


where o22 is the microscopic absorption cross section at 2200 m/sec. The g and s 
factors depend on the shape of the absorption cross section as a function of neutron 
energy. Specifically for nuclides obeying the 1/v law over the entire energy range, 
g = lands = 0. 

The accuracy of the results obtained depends on the input data used which, in 
general, has been taken from the 1958 revision of the Brookhaven Neutron Cross 
Section Compilation. Tables of o and g and s factors in 20 centigrade degree steps 
from 20°C to 760°C are listed for r = 0.03 or r = 0.07. These values of r correspond 
to average parameters appropriate for the moderator and fuel rods respectively of 
the NRX Reactor. In some instances o for r = 0 are given. Elements which follow 
the 1/v law fairly closely in the thermal region do not usually have g values listed. 
The applicability of the compilation is limited to well moderated reactors and to 
thin samples in which self shielding has been neglected. 


K. SHuRE 
Westinghouse Electric Corporation 


Bettis Atomic Power Division 
Pittsburgh, Pennsylvania 


20(T].—Atvin Guiassner, The Thermochemical Properties of the Oxides, Fluorides, 
and Chlorides to 2,500° K, Argonne National Laboratory ANL-5750, U. 8S. 
Government Printing Office, Washington, D. C., 1957, vi + 74 p., 27 cm. Price 
$.45. 


This table gives empirical equations for the thermodynamic properties per mole 
(heat capacity Cp, enthalpy H, entropy S, and gives free energy F) in the form of 
power series in the absolute temperature 7. The heat capacity at constant pressure 
is fitted to the equation: 


d X 10° 

T? 
Only three parameters are evaluated; either c or d is set equal to zero. Integration 
of the heat capacity to give H, S, and F requires two additional constants of inte- 
gration, A and B. The coefficients are tabulated for solid, liquid, and gaseous states 
of the elements (Table I—3 p.), the oxides (Table II—4 p.), the fluorides (Table 
III—5 p.) and the chlorides (Table IV—5 p.). Each of these tables includes, in 
addition, the heat and entropy of the phase transitions, the entropy at 298°K, and 
appropriate references to the source material. Tables V, VI, and VII give the en- 
thalpy and free energies of formation of each substance from the elements at 
298°K, as well as the coefficients Aa, Ab, Ac, Ad, AA, and AB needed to calculate 
values at other temperatures. 

The publication concludes with 21 pages of graphs showing the temperature 
dependence of AF; from 300° to 2500°K. 


C,=a+(bX10°)T + (eX 10°)T + 





Rosert L. Scorr 
University of California 
Los Angeles, California 
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21[W].—Kennetu J. Arrow, Samuet Karun & Hersert Scarr, Studies in the 
Mathematical Theory of Inventory and Production, Stanford University Press, 
Stanford, California, 1958, x + 340 p., 25 em. Price $8.75. 


This work is the initial volume in the Stanford Mathematical Studies in the 
Social Sciences. The reviewer joins the principal authors in recommending careful 
attention to the first two chapters: in Chapter I, Arrow presents a remarkably con- 
cise and enlightening discussion which, more constructively than anything else the 
reviewer has read, relates inventory theory to economics; in Chapter II, this useful 
survey is continued as the principal authors treat common features of many in- 
ventory models after placing them within a realistic framework for decision models. 
The Introduction ends with summaries of results of the remaining three parts: 
Optimal Policies in Deterministic Inventory Processes, Optimal Policies in Stochas- 
tic Inventory Processes, and Operating Characteristics of Inventory Policies. This 
book is judged to devote reasonable attention to computing problems both for cal- 
culation of solutions and for illumination. Reading of individual chapters has con- 
vinced the reviewer that the general promises on computing made by the authors 
on pages 16-19 were honestly kept. The frequent graphs and tables are uniformly 
helpful and pleasing. A bibliography of four pages covering mainly the years 1955- 
1957 is also included. It is to be hoped that subsequent volumes of this Stanford 
Series will push forward into the wide reaches of inventory problems including, for 
example, areas of demand prediction, measures of utility for satisfying differing 
demand patterns, and even seemingly prosaic questions such as how to maintain 
records of extensive inventory systems. In summary, this book is a substantial con- 
tribution to the mathematics of inventory and production problems. Since it pro- 
vides a sound exposition over quite a broad range, it should serve as a valuable 
source for further research. 

W. H. Mariow 


The George Washington University 
Washington, District of Columbia 


22[W, Z).—DanieL D. McCracken, Harotp Weiss, & Tsai-Hwa Lees, Pro- 
gramming Business Computers, John Wiley & Sons, Inc., New York, 1959, 
xvii + 510 p., 24 em. Price $10.25. 


Here is a well-written book about the elements of programming high-speed 
electronic computers. It is particularly written for those people who are interested 
in the programming of management data problems. The book touches upon down- 
to-earth details such as verifying the program accuracy, input and output pro- 
gramming, and rerun techniques. It discusses the advantages and disadvantages of 
machine-aided coding. In general, Programming Business Computers is a compre- 
hensive survey of programming with special emphasis on business data processing. 


ARTHUR SHAPIRO 
Applied Mathematics Laboratory 
David Taylor Model Basin 
Washington, District of Columbia 


23[X].—J. N. Gooprer & P. G. Hones, Jr., Elasticity and Plasticity, v. 1, Surveys 
in Applied Math. Series, John Wiley & Sons, Inc., New York, 1958, ix + 152 p., 
23 em. Price $6.25. 
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The articles of this series are said to be “aimed at a broad, mathematically 
literate audience looking for an up-to-date account of modern progress in applied 
mathematics and an appraisal of future promising research directions.” The first 
author of Volume 1, J. N. Goodier, who contributed the article entitled, ‘“The 
Mathematical Theory of Elasticity,” assumed that the reader was well versed in 
the theory of elasticity, and made no effort to make his article self-contained. He 
contented himself with a brief survey of “those significant recent developments be- 
lieved least known to readers whose first language is English.”” However, even this 
intention is not fully carried out and a list of topics omitted is given at the end of 
the article. Three pages of bibliography are given, which include only those books 
and papers actually discussed or cited. 

The major portion of the discussion deals with work of Russian authors. Great 
stress is given to the work of Muskhelishvili and to investigations inspired by it. 

There is no mention in this article of the application of numerical methods to 
problems in elasticity, aside from a reference to the survey of numerical methods 
in conformal mapping given by G. Birkhoff, D. M. Young and H. Zarantonello, 
in Proc. Symp. Appl. Math, v. 4, 1953, p. 117. 

The second article written by Phillip G. Hodge and entitled, ““The Mathematica! 
Theory of Plasticity,” is practically self-contained and satisfies, very well, the needs 
of the member of the audience described in the first paragraph of this review. 

Chapters 1 to 4 of this article are theoretical in nature, and Chapters 5 to 7 are 
concerned with applications of the theory. In particular, Chapter 5 is a well written, 
moderately exhaustive treatment of the behavior of a simply supported circular 
plate under a uniform normal pressure. In Chapter 6 other problems are discussed 
more briefly in an attempt to illustrate the current state of development in plasticity 
problems. 

Particular attention is paid to significant Russian contributions. Chapter 7 
contains transcription of parts of a report by W. Prager on Russian contributions 
up to 1949, and a section by the author entitled, “Contributions from 1949 to 1955.” 

There is no mention made in this article of the applications of numerical methods 
to problems in plasticity. 

ee 


24[X].—L. V. Kantrorovicn & V. I. Krytov, Approximate Methods of Higher 
Analysis, Translated from the third Russian edition by Curtis D. Benster, Inter- 
science Publishers, Inc., New York, 1958, xv + 681 p., 24 em. Price $17.00. 


In the April 30, 1959 issue of Le Monde, on page 5, there is a description of the 
organization of scientific activities in Russia, in the course of which the following 
remark is made: ‘“‘Contrairement aux Américains, les Russes paraissent parfaite- 
ment au courant de la littérature mondiale.”’ The author is Maurice Letort, ‘“prési- 
dent du comité consultatif de la recherche scientifique et technique.” 

One would like to be indignant, but unfortunately the gibe is deserved. In fact, 
many Americans who visit Russia, or otherwise make contacts with Russian scien- 
tists, are amazed at how up-to-date their acquaintance is with American literature, 
which implies that their own acquaintance with Russian literature is much less so. 
However, the article in Le Monde also provides a partial explanation by describing 
the extensive Russian facilities for translating and abstracting (2000 full time 
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“eollaborators,” plus many part time). It is hardly necessary to comment on the 
meagerness of our own facilities. 

Until recently, even reviews and abstracts of Russian literature were pitifully 
sparse, although here there has been vast improvement. The original Kantorovich- 
Krylov has been known and appreciated by a few Americans, probably largely due 
to informal publicity given it by George Forsythe, who encouraged the making 
and publishing of the present translation. But Mathematical Reviews has no review 
of the second edition, published in 1941, and for the third edition it listed chapter 
headings and remarked, in an unsigned article, only that the edition differed very 
little from the previous one. 

At any rate, we can be grateful to translator and publisher for the present 
volume. The book itself is concerned mainly with the numerical solution of partial 
differential equations, as the title to the first edition (1936) indicated. The first 
chapter deals with expansion in series, both orthogonal and nonorthogonal, with a 
section on the improvement of convergence. Next come methods of solution of 
Fredholm integral equations with applications to the Dirichlet problem. Then comes 
a chapter on difference methods, and one on variational methods. This accounts for 
slightly more than half of the book. There follow two chapters, for a total of nearly 
250 pages, on conformal methods, and finally about 50 pages on Schwarz’s method. 
Throughout, the presentation is extremely readable, with the inclusion of numerous 
examples, but no exercises. Unfortunately there is no index, either, although the 
table of contents is fairly detailed (5 pages). 

In organization the translation deviates from the original only in collecting the 
references at the end, with footnotes referring to author and number. This I con- 
sider to be desirable. In detail the translation is faithful and quite clear. At times 
the phraseology is too faithful for elegance, and on rare occasions the translator is 
even led astray. One such example occurs on page 7: “Just as there, we may separate 
the problem into two, and moreover in each case the conditions are null on two 
sides.”’ While the reader should understand what is meant, there are two faults to 
find here. First, “‘priéem’’ should be translated as “where,” not “and moreover.” 
Second, a condition cannot be null. I confess, I do not understand the construction 
in the original, which is ‘‘usloviya nulevye,” and perhaps the translator can be for- 
given for assuming the adjective to be in predicate form in spite of the ending. 
Perhaps the authors themselves were careless. 

But one could always find fault with details, whereas the important thing is that 
the book is now available to readers of English. Again our thanks to publishers and 
translator. 

A. S. HousEHOLDER 


Oak Ridge National Laboratory 
Oak Ridge, Tennessee 


25(Z|.—Franz L. Aut, Electronic Digital Computers, Academic Press, New York, 
1958, x + 336 p., 23 em. Price $10.00. 


In the preface to his book, Alt addresses himself primarily to “physicists, chem- 
ists, engineers and others in similar occupations who have occasion to require the 
solution of computational problems by means of digital computing machines.” 
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Alt’s primary motivation with respect to this audience is to provide sufficient intro- 
ductory information to improve communication between the “originator of a com- 
puting problem and the team around the machine.” 

The Introduction, Part I, leads the reader gracefully toward Alt’s objectives 
with a well formulated statement of the stages through which a problem goes on 
its way to the number factory. Having outlined the required process of problem 
formulation to problem analysis to programming to coding to machine computa- 
tion, the author proceeds to discuss the process in reverse. Thus, Part II provides 
a functional survey of automatic digital computers and Part III discusses program- 
ming and coding. By the time the object audience reaches more familiar ground in 
Part IV (Problem Analysis) and the statement of computer applications in Part V, 
it has gained the hindsight necessary to make modern computing methods more 
palatable. 

The historical survey of automatic digital computers in Part II is somewhat 
weakened by the selection of the memory type as a principal classification charac- 
teristic. The more significant factor of operation time, or conceptual factors such as 
the stored program are therefore undermined. However, the bulk of Part II ade- 
quately introduces the key factors differentiating one machine from another, that 
is, number representation and memory, arithmetic, control and input-output organs. 
The reviewer felt a lack of graphical presentation of the material in Part II—the 
inclusion of so much data in the same form as the normal prose is likely to lose the 
reader for whom this book is written. 

Part III remains consistent in its reverse discussion of programming and coding, 
by describing coding first and then programming. A 4-address instruction set is de- 
fined and the coding of a simple arithmetic expression and trigonometric function 
is used as a vehicle for explaining coding operations. Single address coding is also 
described in the terminology of IBM manuals. The’ sections on programming 
demonstrate the use of flow charts and define factors and terminology significant, to 
computing tactics. 

Part IV, covering problem analysis, is by far the strongest part of the book. 
It collects, in very readable form, methods of computation used in numerical solu- 
tion of ordinary differential, partial differential, and algebraic equations. As men- 
tioned in the preface, there is no attempt at rigor in this presentation. The discus- 
sion of numerical methods is well seasoned with qualitative comments reflecting 
computational experience and with references to the 210 papers listed in the ap- 
pended bibliography. 

This book cannot by itself serve as a text for the classroom. This book will not 
serve as a reference textbook in the sense of its cataloging completeness. However, 
it is this reviewer’s feeling that Electronic Digital Computers will be found on ref- 
erence shelves for many years, by virtue of its very readable presentation of Alt’s 
extensive experience in high-speed computation. The goal, established in the preface, 
is very adequately achieved. 


GERALD Estrin 
University of California 
Los Angeles, California 
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26|Z|.—_James T. CuLBERTSON, Mathematics and Logic for Digital Devices, D. van 
Nostrand Co., New Jersey, 1958, x + 224 p., 23 em. Price $4.85. 


This text is intended for students with a background of college algebra who plan 
to enter the computer field. As the author points out in his first sentences, it “pre- 
sents neither a course in programming computers nor a mathematical analysis of 
computing mechanisms. Preliminary to these things, it provides a course in some 
mathematics which the student will find useful later”. Specifically, the material 
covered includes some combinatorics and probability, Boolean algebra and proposi- 
tional calculus, with applications to switching circuits. These topics are covered by 
a text written in a highly readable colloquial style, plentifully illustrated by ex- 
amples and by the tremendous total of 1262 numbered exercises. (The latter are 
“dressed-up”, provided with continuity and some attempt at humor in a manner 
which the reviewer finds slightly repellent—but this is a matter of taste. The device 
does have the merit of permitting problems which are essentially identical to appear 
in radically different formulations.) 

The unifying concept of the first part of the text is that of the neuron model. 
Various successive refinements of the receptor-central-effector system are intro- 
duced, leading to adders and other complicated input-output systems. Chapter I is 
introductory, presenting the summation and product notations and the ideas of 
an algorithm and iterative approximation. Unfortunately, two of the five examples 
of the summation notation are incorrect, and there is no clear distinction between 
an algorithm and an iteration. After correctly defining the former as concluding 
in a finite number of steps, the author introduces ““Newton’s algorithm” for the 
square root, which is an example of the iteration process. Further, he states the 
completely false result that B = 4(A + N/A) is always a better approximation 
to ~/N than A. These points will illustrate that statements in the text must be care- 
fully watched; accuracy has frequently been sacrificed to simplicity of statement. 

Chapters II and III present the basic facts about combinations, permutations, 
and elementary frequency probability. Chapter IV is an excellent elementary ac- 
count of arithmetic in various radix systems and of conversion from one system to 
another. The various mechanical procedures, such as subtraction by complementa- 
tion and division by subtraction, are considered in detail. Except for minor matters 
of choice of language, this chapter may be highly recommended. 

The second portion of the book deals with logic. The larger part of Chapter V 
is an exposition of the syllogistic logic in its full mediaeval pattern, including even 
the vowel notation of the scholastics, and omitting only the traditional mnemonics, 
Barbara, Darii, etc. The reviewer finds this portion of the book utterly astounding. 
It is as if one were to come across a long commentary on De Rarum Natura in a 
text on modern atomic physics. Neither mathematicians nor computer engineers 
use syllogisms; what purpose can this chapter serve? The latter portion discusses 
relations, omitting reflexivity and giving much more stringent definitions of one- 
many and many-one than customary. 

The final three chapters develop successively Boolean algebra, the propositional 
calculus, and the model of the latter in terms of switching circuits. The reduction to 
normal form and inivial simplification are well presented. Too much reliance is 
placed on checking by means of Venn diagrams. The student is not adequately 
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warned that, while a Venn diagram illustration of Boolean inequalities is always 
available, equality in a particular Venn diagram does not necessarily imply general 
equality. 

In summary, the text covers matters of algebra, arithmetic, and logic that stu- 
dents should know before taking advanced courses in logical programming or com- 
ponent design. It may be recommended as a text or for collateral reading if the in- 
structor will warn the student of possible pitfalls and inaccuracies. 


J. D. Swirr 
University of California 
Los Angeles, California 


27(Z|.—Depr. or THE Army, Catalog of Commercially Available Automatic Data 
Processing Systems, Department of the Army Pamphlet No. 1-250-4, 1958, 
107 p., 26 em. 


This pamphlet is an excellent compilation and description of automatic data 
processing systems which are commercially available as of July, 1959. It contains 
the descriptions of twenty-five digital data processing systems ranging in size from 
the Bendix G15D and the Royal McBee LGP30 to the UNIVAC 1105 and the 
DATAmatic 1000. A photograph is included with each computer description. 

This compilation is subdivided into categories entitled, respectively, General 
Description, System Components, Programming, Personnel Requirements, and 
Site Preparation. 

The systems component category describes a typical computer configuration 
consisting of central processor, arithmetic unit, input-output control, high speed 
memory, magnetic tape units, paper tape units, card readers and punches, and high 
speed printers. This category lists some pertinent characteristics such as word length, 
numeric characters per word, timing, pulse repetition rate, size of memory, check- 
ing, and error correcting features. A rather complete list of specifications and charac- 
teristics of input and output media is included. The instruction word structure is 
also presented. 

The personnel requirements category recommends a programming and operating 
complement of personnel but excludes maintenance personnel. Manufacturer’s 
training of operators and programmers is discussed, with an option of training at the 
manufacturer’s premises or at the installation site. 

The over-all floor space, floor loading, and air conditioning requirements are 
also given in the site preparation category. 

However, the major contributions of this pamphlet are the tables of cost, power 
requirements, and physical characteristics of each unit. These tables also present 
rental, purchase, and maintenance costs. 

This reviewer considers this pamphlet an excellent guide to all computer users 
who require a ready reference in the field. 


ALEXANDER C. ROSENBERG 
Applied Mathematics Laboratory 
David Taylor Model Basin 
Washington, District of Columbia 
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28(Z|.—Martin GARDNER, Logic Machines and Diagrams, McGraw-Hill Book 
Company, Inc., New York, 1958, ix + 157 p., 23 em. Price $5.00. 


Only on rare occasions does one come across a book written on a technical sub- 
ject which is entertaining as well as informative. It is also unusual to find a book on 
the mechanistic aspects of formal logic which does not begin with either Venn 
diagrams or a description of Boolean algebra. The author has presented an historical 
survey of the subject in a somewhat narrative fashion, beginning with an almost 
complete biography of Ramon Lull and ending with speculations on the future of 
logic machines. The ‘‘References” after each chapter are considerably more than 
just references, and make as interesting reading as the text. The book is by no means 
devoid of the author’s opinions and no attempt has been made at concealment. In 
fact, it is amusing to note that, even though some of the artifacts and methods de- 
scribed in the book are treated somewhat modestly by the author, he cannot resist 
the temptation to devote a chapter to a method which he, himself, has devised. 
This, I am sure, is understandable to any person who has worked in the field. 

Persons interested in the field of logic, either as a subdivision of philosophy or 
as an aid to digital computer design, will find the book well worth its reading time. 
My only complaint is that he has not included the work done in the area of com- 
puter design, which could be interpreted as legitimate subject matter under this 
title. 

NeEtson T. GRISAMORE 
The George Washington University 
Washington, District of Columbia 








TABLE ERRATUM 


273.—E. JAHNKE & F. Emp, Tables of Functions with Formulae and Curves. 4th 
edition, Dover Publications, New York, 1945. 


On p. 99, the caption for Fig. 55 includes values of w, w’, k, k’, e1 , e2 , es , go and 
gs associated with the Weierstrassian @ function. It is, of course, impossible to tell 
with certainty which of these numbers were intended to be assumed as given 
exactly. A reasonable inference, however, is that k = 0.8 and e, — e; = 1 were 
intended as given. If this be the case, the values of e; and g; are in error. The 3D 
values of e; and g; should, then, be —0.547 and —0.093 respectively (indeed, 
es; = —0.5467 and g; = —0.09252, to 4S). 

Tuomas H. SouTHARD 


University of California 
Los Angeles, California 


EpitoriaL Note: Under the assumptions stated in the preceding communication, the 
values of w and w’ are, respectively, 1.995305 and 1.750754: to 6D, which when rounded to two 
decimal places appear to be exact numbers, as one might erroneously infer from the values 
shown to that accuracy in the caption under discussion. 


J. W. W. 
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CORRIGENDA 


Danie. C. Frevper, “A note on summation formulas of powers of roots,’’ MT'AC, v. 12; 
1958, p. 197. 

The term 27a,a2*a;/ao* has been omitted from the expression for Ss. The fourth line in this 
expression should accordingly be corrected to read 


+ (9a;7ag + 27a;2a2a; + 27a;2a20, + 27a,a2a;? + Zayar*a, + Ya2*az)/ao*. 


G. P. M. Hese.pen 


Her Masesty’s Nautica, ALmanac Orrice, Subtabulation, A Companion Booklet to Inter- 
polation and Allied Tables, MTAC, Review 18, v. 13, 1959, p. 127-129. 


for read 
p. 128, first displayed formula pd} pd 
p. 128, last displayed formula ber tr 
B;(rv)C B;(rv)c 
ve + yi ve! + v1‘ 
p. 129, first displayed formula 3 5 


Joun Topp 
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